diff --git a/main.cpp b/main.cpp
index 923dce5..9cb7f94 100644
--- a/main.cpp
+++ b/main.cpp
@@ -1,208 +1,210 @@
 #include <iostream>
 #include <fstream>
 #include <Eigen/Dense>
 #include <vector>
 #include <cmath>
 #include "ODElibrary/ForwardEuler.h"
 #include "ODElibrary/RungeKutta.h"
 #include "ODElibrary/AdamsBashforth.h"
 #include "ODElibrary/AdamsMoulton.h"
 #include "ODElibrary/BDF.h"
 
 
 
 /**
  * This function will correspond to the rhs function g(t), where:
  * @param t: time at which it has to be computed.
  * It is designed to compute the values at a single specific time.
  */
 Eigen::VectorXd g_function(double t);
 
 /**\brief This function corresponds to another possible rhs function g(t)  (more details below).
  */
 Eigen::VectorXd ForcedOscillations(double t);
 
 
 int main(int argc, char *argv[]){
 
     std::cout << "Number of command line arguments = "
               << argc << "\n";
     for (int i=0; i<argc; i++) {
         std::cout << "Argument " << i << " = " << argv[i];
         std::cout << "\n";
     }
 
 
     Eigen::VectorXd initialValues;
     double parameters[11];
     std::fstream inputfile;
     inputfile.open(argv[1]);
     std:: string s;
     for (int i=0; i<11;i++) {
         inputfile >> s >> parameters[i];
         if (i==7){
             int dimension=(int)(parameters[6]);
             initialValues.resize(dimension);
             initialValues[0]=parameters[7];
             int j=1;
             while(j<dimension){
                 inputfile>>initialValues[j];
                 j++;
             }
         }
 
     }
 
     inputfile.close();
 
 
 
     std::vector<AbstractSolver*> solvers;
 
     int solver_name=(int)(parameters[0]);
 
     switch(solver_name){
         case 1: {         //Euler
             auto solver = new ForwardEuler();
             solvers.push_back(solver);
             break;
         }
         case 2: {         //RK
             auto solver = new RungeKutta();
             solver->setMethodOrder((int)(parameters[10]));
             solvers.push_back(solver);
             break;
         }
         case 3: {         //Adams Bashforth
             auto solver = new AdamsBashforth();
             solver->setStepsAlgo((int)(parameters[9]));
             solvers.push_back(solver);
             break;
         }
         case 4: {          // Adams Moulton
             auto solver = new AdamsMoulton();
             solver->setStepsAlgo((int)(parameters[9]));
             solvers.push_back(solver);
             break;
         }
         case 5: {        //BDF
             auto solver = new BDF();
             solver->setStepsAlgo((int)(parameters[9]));
             solvers.push_back(solver);
             break;
         }
         default:
             std::cerr<<"Method not known. Set a number from 1 to 5.\n";
 
     }
 
     AbstractSolver* solver=solvers[0];
     solver->setInitialTime(parameters[1]);
     solver->setFinalTime(parameters[2]);
     int discretizationMethod=parameters[5];
     if(discretizationMethod==0)
         solver->setNumberOfSteps(parameters[3]);
     else
         solver->setStepSize(parameters[4]);
 
     solver->setInitialValue(initialValues);
 
-
     Eigen::VectorXd (*p_function)(double t);
-    p_function=&g_function;
     Eigen::MatrixXd m(2,2);
-    m(0,0)=1.0, m(0,1)=0.0,m(1,0)=0.0, m(1,1)=-2.0;
-    double springStiffness=-0.001;
-    Eigen::MatrixXd p(2,2); // from spring equation
-    p(0,0)=0, p(0,1)=1,p(1,0)=-springStiffness,p(1,1)=0.03;
-    Eigen::VectorXd (*p_forcedoscillations)(double t);
-    p_forcedoscillations=&ForcedOscillations;
+    if((int)(parameters[8])==1) {
+        p_function = &g_function;
+        m(0,0)=1.0, m(0,1)=0.0,m(1,0)=0.0, m(1,1)=-2.0;
+    }
+    else{
+        p_function=&ForcedOscillations;
+        double springStiffness=-0.001;
+        m(0,0)=0, m(0,1)=1,m(1,0)=-springStiffness,m(1,1)=0.3;
+    }
+
 
-    solver->setRightHandSide(p,p_forcedoscillations);
+    solver->setRightHandSide(m,p_function);
 
     std::ofstream outputfile;
     outputfile.open("solution.dat");
     solver->solve(outputfile);
     outputfile.close();
 
 
 
 /*
     Eigen::VectorXd y0(2);
     y0[0]=1.,y0[1]=1.;
     /**
      * Another possible problem that can be used to test this library describes the behaviour of a spring with a mass attached to it,
      * with an additional term representing forced oscillations. The matrix describing it is here defined, along with the
      * initial vector used to solve it, y_1(0)=0, y_2(0)=0. At t=1, the exact solutions of the vectorial differential equation
      * is y_1(1)=0.182532458, y_2(1)=0.046017418.
 
     double springStiffness=0.1;
     Eigen::MatrixXd p(2,2); // from spring equation
     p(0,0)=0, p(0,1)=1,p(1,0)=-springStiffness,p(1,1)=0.3;
     Eigen::VectorXd (*p_forcedoscillations)(double t);
     p_forcedoscillations=&ForcedOscillations;
     Eigen::VectorXd yStart(2);
     yStart[0]=0.,yStart[1]=0.;
 
 
 
     /**
      * In the following part of the code, different solvers can be called to solve the ODE problem at hand.
      * By leaving uncommented -ONLY- the line corresponding to the method which one wants to apply, the problem
      * is solved with it. Parameters can be played with.
 
     double initialTime=0;
     double finalTime=1;
     int nSteps=100;
     int nStepsMultistep=5;
     int orderRK=2;
 
     /**\brief Change y0 to yStart to integrate the spring problem instead of the exponential test problem.
 
     //ForwardEuler solverq=ForwardEuler(initialTime,finalTime,nSteps,yStart);
     //RungeKutta solverq=RungeKutta(initialTime,finalTime,nSteps,orderRK,yStart);
     //AdamsBashforth solverq=AdamsBashforth(initialTime,finalTime,nSteps,nStepsMultistep,yStart);
     //AdamsMoulton solverq=AdamsMoulton(initialTime,finalTime,nSteps,nStepsMultistep,yStart);
     //BDF solverq=BDF(initialTime,finalTime,nSteps,nStepsMultistep,yStart);
 
     /**
      * \brief The rhs function is set to the desired value. Any matrix and any pointer to a function returning a
      * Eigen::VectorXd can be passed as parameters. The solver will then use them to integrate the differential equations.
      * Setting solver.setRightHandSide(p,p_forcedoscillations), the spring problem described above can be solved.
 
     //solverq.setRightHandSide(m,p_function);
     //solverq.setRightHandSide(p,p_forcedoscillations);
 
     /**
      * The problem is solved and the solution is displayed on std::cout, since this is the default of the 'solve' method.
      * By calling solver.solve(solutionfile) the computed values will instead be written to the file 'solutionfile'.
 
     std::ofstream solutionfile;
     solutionfile.open("AdamsBashforth2_solution.txt");
     //solverq.solve();
     solutionfile.close();
 
 */
 }
 
 /**
  * Implementation of the rhs functions g. The entries of the result are set to zero so that the accuracy of the simple integration
  * performed in this main can be easily checked. It can be changed to whatever is desired by the user. The return value
  * does not have to be two-dimensional, the function can be modified to return a vector of whatever length.
  * It is set to two dimensions just for demonstration purposes.
  * @param t : time at which the function is computed
  * @return  : An Eigen::VectorXd storing the value of g(t).
  */
 Eigen::VectorXd g_function(double t){
     Eigen::VectorXd result(2);
     result[0]=0;
     result[1]=0;
     return result;
 }
 
 Eigen::VectorXd ForcedOscillations(double t){
     Eigen::VectorXd result(2);
     result[0] = 0;
     result[1] = cos(2*M_PI*t); // derivative of velocity
     return result;
 };
\ No newline at end of file