diff --git a/ODElibrary/AbstractSolver.cpp b/ODElibrary/AbstractSolver.cpp
index 72e7a90..998d524 100644
--- a/ODElibrary/AbstractSolver.cpp
+++ b/ODElibrary/AbstractSolver.cpp
@@ -1,148 +1,156 @@
 //
 // Created by lionel on 25.11.19.
 //
 
 #include <cassert>
 #include <sstream>
 #include <iostream>
 #include "AbstractSolver.h"
 
 
 
 AbstractSolver::AbstractSolver(){
     stepSize=0;
     initialTime=0;
     finalTime=0;
     dimension=0;
     nSteps=0;
 };
 
+AbstractSolver::AbstractSolver(double t0,double tf,int n,Eigen::VectorXd& y0){
+    initialTime = t0;
+    finalTime = tf;
+    nSteps = n;
+    stepSize = (tf - t0) / n;
+    dimension=y0.size();
+    initialValue = y0;
+}
 
 AbstractSolver::~AbstractSolver(){};
 
 void AbstractSolver::setInitialTime(const double t0) {
     initialTime = t0;
 }
 
 void AbstractSolver::setFinalTime(const double tf) {
     finalTime = tf;
 }
 
 void AbstractSolver::setStepSize(const double h) {
     checkStepSize(h);
     if (stepSize!=0){ // need to warn user that the stepsize was changed to avoid confusion
         std::cerr << "[WARNING] : Step size " << stepSize << " was already assigned, new step size is now " << h << std::endl;
     }
     stepSize = h;
 }
 
 void AbstractSolver::checkStepSize(double h) {
     if (h <= 0) { // negative step size will make our solver fail
         throw std::domain_error("Step size h must be strictly positive");
     }
 }
 
 
 void AbstractSolver::checkTimes() {
     if (finalTime <= initialTime){ // final time must be larger for the solve method to work
         throw std::domain_error("Initial time is not smaller than the final time.");
     }
 }
 
 void AbstractSolver::setNumberOfSteps(int n){
     if (n < 1){
         throw std::domain_error("Number of step must be strictly positive ");
     } else {
         double h = (finalTime - initialTime)/ (double) n; // linearly spaced steps
         setStepSize(h);
     }
     nSteps=n;
 }
 
 void AbstractSolver::checkSpaceDimension(int dim) {
     if (dim < 1) {
         throw std::domain_error("Dimension must be greater than 0");
     } else{
         if (dimension == 0) { // the dimension was never specified if it is 0, so we set it to dim
             dimension = dim;
         } else if (dimension != dim){ // return an error with the conflicting values of dimension
             std::stringstream errorMessage;
             errorMessage << "Invalid dimension : " << dim << ", it was already assigned " << dimension << ".";
             throw std::length_error(errorMessage.str());
         } else { // do nothing if the values of dimension match
             ;
         }
     }
 }
 void AbstractSolver::setInitialValue(Eigen::VectorXd & y0) {
     // check it is a valid input
     checkSpaceDimension(y0.size());
     initialValue = y0;
 }
 
 void AbstractSolver::setRightHandSide(Eigen::MatrixXd & rhsMatrix, Eigen::VectorXd (*g_)(double t)){
     // check the right hand side matrix is valid before storing it
     try {
         checkSpaceDimension(rhsMatrix.rows());
         checkSpaceDimension(rhsMatrix.cols());
     } catch (const std::length_error& e){
         std::stringstream errorMessage; // creates a more complex message to help the user correcting his error
         errorMessage << "The linear Matrix has dimension ("<< rhsMatrix.rows() <<","<< rhsMatrix.cols() <<
                                                         ") which is invalid : "<< e.what();
         throw std::length_error(errorMessage.str());
     }
     A = rhsMatrix;
 
     // check if the function is callable and that it returns object of desired dimension
     try {
         // checks if no exception happens for final and initial time, hope no exception will happen for in between times
         auto testVector = g_(initialTime);
         auto testVectorCopy = g_(finalTime);
     } catch(const std::exception& e) {
         std::stringstream errorMessage;
         errorMessage << "Function g(t) could not be evaluated, an exception was caught with message : " << e.what();
         throw std::runtime_error(errorMessage.str());
     } catch (...) { // any other type of error will make our solver fail
         throw std::runtime_error("Function g(t) could not be evaluated. No exception was caught.");
     }
     // checks the output of g is in the required format
     try {
         checkSpaceDimension(g_(initialTime).size());
     } catch (const std::length_error& e){
         std::stringstream errorMessage;
         errorMessage << "Function g(t) has invalid dimension in output : " << e.what();
         throw std::length_error(errorMessage.str());
     }
     g = g_;
 }
 
 double AbstractSolver::getInitialTime() const{
     return initialTime;
 }
 
 double AbstractSolver::getFinalTime() const{
     return finalTime;
 }
 
 double AbstractSolver::getStepSize() const{
     return stepSize;
 }
 
 int AbstractSolver::getSpaceDimension() const {
     return dimension;
 }
 
 void AbstractSolver::solve(std::ostream &stream) {
     auto y = initialValue;
     double t = initialTime;
     assert(initialTime < finalTime && stepSize > 0); // ensures next while loop will not be infinite
     stream << t << " " << y.transpose() << std::endl; //write the initial time
     while(t< finalTime){
         step(y,t); // step updates the y values
         t += stepSize;
         stream << t << " " << y.transpose() << std::endl;
     }
 }
 
 
 
diff --git a/ODElibrary/AbstractSolver.h b/ODElibrary/AbstractSolver.h
index 4297df0..668c1c1 100644
--- a/ODElibrary/AbstractSolver.h
+++ b/ODElibrary/AbstractSolver.h
@@ -1,104 +1,105 @@
 //
 // Created by lionel on 25.11.19.
 //
 
 #ifndef ODELIBRARY_ABSTRACTSOLVER_H
 #define ODELIBRARY_ABSTRACTSOLVER_H
 
 
 #include <vector>
 #include <Eigen/Dense>
 
 
 
 /** This is an abstract class for the solver of differential equations.
  * It provides a lot of method for setting up a solver.
  * The only pure virtual method that needs to be overriden by the daughter classes is the step method.
  */
 class AbstractSolver {
 protected:
     // variables for the setting
     double stepSize;
     double initialTime;
     double finalTime;
     // dimension of the differential system (size of y)
     int dimension;
     int nSteps;
     Eigen::VectorXd initialValue; // y_0
     /// function that throws exceptions if the dim is not the one of the solver
     void checkSpaceDimension(int dim);
     /// Checks that the time values are different and t0 < tf
     void checkTimes();
     /// Checks that the value of the step size is allowed for our algorithm
     void checkStepSize(double h);
 
     // variables for the rhs of the equation
     Eigen::MatrixXd A; // Linear coefficients of the differential system
     Eigen::VectorXd (*g)(double t); // function g depending on time
 
     /**
      * Step function will perform a step at a given time with a given value associated to that time.
      * It must update the y value in place.
      * It is purely virtual because it is specific to every algorithm.
      * @param y : Eigen vector, the current value of the equation
      * @param t : double, the current time where the step must be performed
      */
     virtual void step(Eigen::VectorXd &y, double t) = 0;
 
 public:
     /// Default constructor, sets everything to 0
     AbstractSolver();
+    AbstractSolver(double t0,double tf,int n,Eigen::VectorXd& y0);
     virtual ~AbstractSolver();
     /// sets the initial value of the time axis
     void setInitialTime(double t0);
     /// sets the final value of the time
     void setFinalTime(double tf);
     /**
      * Sets the starting value at t0.
      * @param y0
      */
     void setInitialValue(Eigen::VectorXd & y0);
     /**
      * Sets the step size, returns an exception if the step size is not strictly positive.
      * If the step size was already assigned it will print a warning in the error channel.
      * @param h : step size
      */
     void setStepSize(double h);
     /**
      * Sets the number of steps that are used for the algorithm
      * @param n : number of step, an integer greater than one
      */
     void setNumberOfSteps(int n);
     /**
      * Sets the right hand size of the equation. The equation is of type y' = A * y + g(t).
      * It will also perform checks that the dimensions correspond with the other values of the problem and check that
      * the function can be evaluated properly.
      * @param rhsMatrix : a Eigen Matrix containing the linear part of the equation
      * @param g_ : the time dependent function, returning an Eigen vector object given the time
      */
     virtual void setRightHandSide(Eigen::MatrixXd & rhsMatrix, Eigen::VectorXd (*g_)(double t));
 
     /// Getter for the initial time value
     double getInitialTime() const;
     /// Getter for the final time value
     double getFinalTime() const;
     /// Getter for the step size value
     double getStepSize() const;
     /// Getter for the dimension value
     int getSpaceDimension() const;
 
     /**
      * The solve method will iterate on the time domain of the function, if it has been specified before.
      * Outputs at every line space separated :
      *      1. the time
      *      2. all elements of the vector y for the given time
      * The output of the solver will be written on the ostream object that is specified in the arguments, default will
      * print to the terminal (std::cout).
      * @param ostream: std::ostream type object on which the output will be printed, default : std::cout
      */
     virtual void solve(std::ostream &stream = std::cout);
 
 };
 
 
 #endif //ODELIBRARY_ABSTRACTSOLVER_H
\ No newline at end of file
diff --git a/ODElibrary/AdamsBashforth.cpp b/ODElibrary/AdamsBashforth.cpp
index 642172e..39a26f8 100644
--- a/ODElibrary/AdamsBashforth.cpp
+++ b/ODElibrary/AdamsBashforth.cpp
@@ -1,55 +1,66 @@
 #include "AdamsBashforth.h"
 //#include "MultistepAbstractSolver.h"
 
 AdamsBashforth::AdamsBashforth():MultistepAbstractSolver(){
 };
 
 AdamsBashforth::AdamsBashforth(double t0,double tf,int n,int steps,Eigen::VectorXd& y0):MultistepAbstractSolver(t0,tf,n,steps,y0){
     convergenceRate=steps;
     /**\brief
      * According to the number of steps n desired, the coefficients defining the algorithm are stored in bCoefficients.
      * For implementation purposes, the order of the coefficients is reversed with respect to the classic formula.
      */
-    bCoefficients.resize(steps);
-    switch(steps){
+     setBCoeff(steps);
+};
+
+void AdamsBashforth::setStepsAlgo(int n){
+    stepsAlgo=n;
+    convergenceRate=n;
+    setBCoeff(n);
+}
+
+void AdamsBashforth::setBCoeff(int n) {
+    bCoefficients.resize(n);
+    switch (n) {
         ///When steps=1, Adams Bashforth multistep algorithm actually coincides with Forward Euler.
         case 1:
-            bCoefficients[0]=1;
+            bCoefficients[0] = 1;
             break;
         case 2:
-            bCoefficients[0]=-0.5,bCoefficients[1]=1.5;
+            bCoefficients[0] = -0.5, bCoefficients[1] = 1.5;
             break;
         case 3:
-            bCoefficients[0]=5./12. ,bCoefficients[1]=-16./12.,bCoefficients[2]=23./12.;
+            bCoefficients[0] = 5. / 12., bCoefficients[1] = -16. / 12., bCoefficients[2] = 23. / 12.;
             break;
         case 4:
-            bCoefficients[0]=-9./24.,bCoefficients[1]=37./24.,bCoefficients[2]=-59./24.,bCoefficients[3]= 55./24.;
+            bCoefficients[0] = -9. / 24., bCoefficients[1] = 37. / 24., bCoefficients[2] =-59. / 24., bCoefficients[3] = 55. / 24.;
             break;
         default:
-            std::cerr<<"The number of steps has to be an integer between 1 and 4."<<"\n";
+            std::cerr << "The number of steps has to be an integer between 1 and 4." << "\n";
     }
-};
+}
+
 
 void AdamsBashforth::step(Eigen::VectorXd &y, double t){
     auto it=rhsPreviousSteps.begin();
     int count=0;
     /**\brief
      * The following cycle, which coincides with one integration step, defines the update of y.
      */
     for (it=rhsPreviousSteps.begin();it!=rhsPreviousSteps.end();it++) {
         y +=(*it)*stepSize*bCoefficients[count]; //reversed coefficients in bCoefficients!
         count++;
     }
     /**\brief
      * Having computed yNext, the rhs function can subsequently be evaluated at tNext.
      */
     Eigen::VectorXd rhsNext = A *y+g(t+stepSize);
     /**\brief
      * Such value is then stored in rhsPreviousStep. The value at the beginning of the list is then discarded,
      * since it is not needed for the next integration step.
      */
     rhsPreviousSteps.push_back(rhsNext);
     rhsPreviousSteps.pop_front();
 }
 
 
diff --git a/ODElibrary/AdamsBashforth.h b/ODElibrary/AdamsBashforth.h
index 940d722..9696253 100644
--- a/ODElibrary/AdamsBashforth.h
+++ b/ODElibrary/AdamsBashforth.h
@@ -1,21 +1,23 @@
 #ifndef ODELIBRARY_ADAMSBASHFORTH_H
 #define ODELIBRARY_ADAMSBASHFORTH_H
 
 
 #include <iostream>
 #include <Eigen/Dense>
 #include <list>
 #include <vector>
 #include "MultistepAbstractSolver.h"
 
 class AdamsBashforth:public MultistepAbstractSolver{
 public:
     AdamsBashforth();
     AdamsBashforth(double t0,double tf,int n,int steps,Eigen::VectorXd& y0);
 
+    void setStepsAlgo(int n);
+    void setBCoeff(int n);
     void step(Eigen::VectorXd &y, double t) override;
 };
 
 
 
 #endif
\ No newline at end of file
diff --git a/ODElibrary/AdamsMoulton.cpp b/ODElibrary/AdamsMoulton.cpp
index f68722c..61bce6d 100644
--- a/ODElibrary/AdamsMoulton.cpp
+++ b/ODElibrary/AdamsMoulton.cpp
@@ -1,82 +1,91 @@
 #include <Eigen/Dense>
 #include "AdamsMoulton.h"
 
 AdamsMoulton::AdamsMoulton():MultistepAbstractSolver(){
 };
 
 AdamsMoulton::AdamsMoulton(double t0,double tf,int n,int steps,Eigen::VectorXd& y0):MultistepAbstractSolver(t0,tf,n,steps,y0){
     convergenceRate=steps+1;
 };
 
+void AdamsMoulton::setStepsAlgo(int n){
+    stepsAlgo=n;
+    convergenceRate=n+1;
+    setBCoeff();
+}
 
-void AdamsMoulton::setRightHandSide(Eigen::MatrixXd & rhsMatrix, Eigen::VectorXd (*g_)(double t)){
-    /**\brief
-     * At first, the function of the abstract mother class AbstractSolver is called. It initializes the actual rhs function.
-     */
-    AbstractSolver::setRightHandSide(rhsMatrix, *g_);
-    /**\brief
-     * This Matrix is needed to build the updateMatrix which defines Adams Moulton method.
-     */
-    Eigen:: MatrixXd I=Eigen::MatrixXd::Identity(dimension,dimension);
-    /**
-     * \brief Depending on the number of steps characterizing the algorithm, the vector bCoefficients is initialiized with the proper
-     * coefficients defining the updates of the method.
-     */
+void AdamsMoulton::setBCoeff() {
     bCoefficients.resize(stepsAlgo+1);
     switch(stepsAlgo){
         case 1:
             bCoefficients[0]=0.5,bCoefficients[1]=0.5;
             break;
         case 2:
             bCoefficients[0]= -1./12.,bCoefficients[1]=8./12.,bCoefficients[2]=5./12.;
             break;
         case 3:
             bCoefficients[0]=1./24.,bCoefficients[1]=-5./24.,bCoefficients[2]=19./24.,bCoefficients[3]= 9./24.;
-        break;
+            break;
         case 4:
             bCoefficients[0]=-19./720.,bCoefficients[1]=106./720.,bCoefficients[2]= -264./720.,bCoefficients[3]=646./720.,bCoefficients[4]=251./720.;
             break;
         default:
             std::cerr<<"The number of steps has to be an integer between 1 and 4."<<"\n";
 
     }
+}
+
+void AdamsMoulton::setRightHandSide(Eigen::MatrixXd & rhsMatrix, Eigen::VectorXd (*g_)(double t)){
+    /**\brief
+     * At first, the function of the abstract mother class AbstractSolver is called. It initializes the actual rhs function.
+     */
+    AbstractSolver::setRightHandSide(rhsMatrix, *g_);
+    /**\brief
+     * This Matrix is needed to build the updateMatrix which defines Adams Moulton method.
+     */
+    Eigen:: MatrixXd I=Eigen::MatrixXd::Identity(dimension,dimension);
+    /**
+     * \brief Depending on the number of steps characterizing the algorithm, the vector bCoefficients is initialiized with the proper
+     * coefficients defining the updates of the method.
+     */
+    setBCoeff();
     /**\brief
      * The matrix involved in the computation of the next step value of y is computed.
      */
     int s=bCoefficients.size();
     double d=bCoefficients[s-1]*stepSize;
     updateMatrix=I-d*A;
 };
 
 void AdamsMoulton::step(Eigen::VectorXd &y, double t){
     /**
      * \brief This vector will correspond to the right hand side of the system to be solved in order to perform one step of the algorithm,
      * updateMatrix*yNext=rhsVector. The vector is then suitably computed.
      */
     Eigen::VectorXd rhsVector(dimension);
     int s=bCoefficients.size();
     rhsVector=y+stepSize*bCoefficients[s-1]*g(t+stepSize);
     int count=0;
     auto it=rhsPreviousSteps.begin();
     for(it=rhsPreviousSteps.begin();it!=rhsPreviousSteps.end();it++){
         rhsVector+=(*it)*bCoefficients[count]*stepSize;
         count++;
     }
     /**
      * \brief The solution at the next time step is computed.
      */
     Eigen::VectorXd yNext=updateMatrix.lu().solve(rhsVector);
     /**
      * \brief The value of the rhs function defining the ODE problem is defined and stored in rhsPreviousSteps.
      */
     Eigen::VectorXd rhsNext=A*yNext+g(t+stepSize);
     rhsPreviousSteps.push_back(rhsNext);
     /**
      * \brief The first value of the list rhsPreviousSteps is not needed for the next integration step, so it is discarded.
      */
     rhsPreviousSteps.pop_front();
     /**
      * \brief The variable y is updated.
      */
     y=yNext;
 };
\ No newline at end of file
diff --git a/ODElibrary/AdamsMoulton.h b/ODElibrary/AdamsMoulton.h
index b6e1d24..7cbded9 100644
--- a/ODElibrary/AdamsMoulton.h
+++ b/ODElibrary/AdamsMoulton.h
@@ -1,48 +1,50 @@
 #ifndef ODELIBRARY_ADAMSMOULTON_H
 #define ODELIBRARY_ADAMSMOULTON_H
 
 #include <iostream>
 #include <Eigen/Dense>
 #include <list>
 #include <vector>
 #include "MultistepAbstractSolver.h"
 
 class AdamsMoulton:public MultistepAbstractSolver{
 private:
     /**
      * When applying an implicit multistep method to solve a linear ODE system, at each iteration the update on the
      * y value is performed solving a suitable linear system. The matrix that defines this system is constant through-out the
      * iterations, whereas the rhs vector of the linear system updateMatrix*yNext=rhsVector is computed at each iteration.
      */
     Eigen::MatrixXd updateMatrix;
 public:
     AdamsMoulton();
     AdamsMoulton(double t0,double tf,int n,int steps,Eigen::VectorXd& y0);
 
+    void setStepsAlgo(int n);
+    void setBCoeff();
     /**
      * This method overrides the setRightHandSide method of the abstract mother class. This operation is needed for implicit multistep
      * methods, because the matrix employed at each iteration for the update of y, updateMatrix, has to be computed and stored
      * as a variable of this class.
      * @param rhsMatrix: Matrix defining the ODE linear system
      * @param g_: function of t, part of the rhs expression.
      */
-    void setRightHandSide(Eigen::MatrixXd & rhsMatrix, Eigen::VectorXd (*g_)(double t));
+    void setRightHandSide(Eigen::MatrixXd & rhsMatrix, Eigen::VectorXd (*g_)(double t)) override;
     /**
      * This function defines the update which is performed at each step of Adams Moulton method.
      * @param y : reference to the last computed value of y
      * @param t : time corresponding to the last computed value of y. The function computes y at time t+stepSize.
      */
     void step(Eigen::VectorXd &y, double t) override;
 };
 
 
 
 
 
 
 
 
 
 
 
 #endif
\ No newline at end of file
diff --git a/ODElibrary/BDF.cpp b/ODElibrary/BDF.cpp
index 7b1ae13..a31abc7 100644
--- a/ODElibrary/BDF.cpp
+++ b/ODElibrary/BDF.cpp
@@ -1,81 +1,91 @@
 #include <Eigen/Dense>
 #include "BDF.h"
 
 BDF::BDF():MultistepAbstractSolver(){
 };
 
 BDF::BDF(double t0,double tf,int n,int steps,Eigen::VectorXd& y0):MultistepAbstractSolver(t0,tf,n,steps,y0){
     convergenceRate=steps;
 };
 
 
-void BDF::setRightHandSide(Eigen::MatrixXd & rhsMatrix, Eigen::VectorXd (*g_)(double t)){
-    /**\brief
-     * At first, the function of the abstract mother class AbstractSolver is called. It initializes the actual rhs function.
-     */
-    AbstractSolver::setRightHandSide(rhsMatrix,*g_);
-    /**\brief
-     * This Matrix is needed to build the updateMatrix which defines Adams Moulton method.
-     */
-    Eigen:: MatrixXd I=Eigen::MatrixXd::Identity(dimension,dimension);
-    /**
-     * \brief Depending on the number of steps characterizing the algorithm, the vector bCoefficients is initialiized with the proper
-     * coefficients defining the updates of the method.
-     */
+void BDF::setStepsAlgo(int n){
+    stepsAlgo=n;
+    convergenceRate=n;
+    setBCoeff();
+}
+
+void BDF::setBCoeff() {
     bCoefficients.resize(stepsAlgo+1);
     switch(stepsAlgo){
         case 1:
             bCoefficients[0]=-1.,bCoefficients[1]=1.;
             break;
         case 2:
             bCoefficients[0]= 0.5,bCoefficients[1]=-2.,bCoefficients[2]=1.5;
             break;
         case 3:
             bCoefficients[0]=-1./3.,bCoefficients[1]=1.5,bCoefficients[2]= -3.,bCoefficients[3]=11./6.;
             break;
         case 4:
             bCoefficients[0]=1./4.,bCoefficients[1]=-4./3.,bCoefficients[2]= 3.,bCoefficients[3]=-4.,bCoefficients[4]= 25./12.;
             break;
         case 5:
             bCoefficients[0]=-1./5.,bCoefficients[1]=5./4.,bCoefficients[2]= -10./3.,bCoefficients[3]=5.,bCoefficients[4]=-5.,bCoefficients[5]=137./60.;
             break;
         default:
             std::cerr<<"The number of steps has to be an integer between 1 and 5."<<"\n";
 
     }
+}
+
+void BDF::setRightHandSide(Eigen::MatrixXd & rhsMatrix, Eigen::VectorXd (*g_)(double t)){
+    /**\brief
+     * At first, the function of the abstract mother class AbstractSolver is called. It initializes the actual rhs function.
+     */
+    AbstractSolver::setRightHandSide(rhsMatrix,*g_);
+    /**\brief
+     * This Matrix is needed to build the updateMatrix which defines Adams Moulton method.
+     */
+    Eigen:: MatrixXd I=Eigen::MatrixXd::Identity(dimension,dimension);
+    /**
+     * \brief Depending on the number of steps characterizing the algorithm, the vector bCoefficients is initialiized with the proper
+     * coefficients defining the updates of the method.
+     */
+    setBCoeff();
     /**\brief
      * The matrix involved in the computation of the next step value of y is computed.
      */
     int s=bCoefficients.size();
     double d=bCoefficients[s-1];
     updateMatrix=d*I-stepSize*A;
 };
 
 void BDF::step(Eigen::VectorXd &y, double t){
     /**
      * \brief This vector will correspond to the right hand side of the system to be solved in order to perform one step of the algorithm,
      * updateMatrix*yNext=rhsVector. This vector is progressively assembled with the suitable expressions.
      */
     Eigen::VectorXd rhsVector(dimension);
     rhsVector=stepSize*g(t+stepSize);
     int count=0;
     /**
      * \brief The vector yInitialConditions originally contains the missing initial conditions computed in order to initialize the algorithm.
      * In the BDF method, this vector is then reused to store the last values of y, needed to perform the next step update.
      */
     auto it=yInitialConditions.begin();
     for(it=yInitialConditions.begin();it!=yInitialConditions.end();it++){
         rhsVector-=(*it)*bCoefficients[count];
         count++;
     }
     /**
      * \brief The solution at the next time step is computed and stored in the list yInitialCondition.
      */
     Eigen::VectorXd yNext=updateMatrix.lu().solve(rhsVector);
     y=yNext;
     yInitialConditions.push_back(yNext);
     /**
      * \brief The first value of the list is not needed for the next integration step, so it is discarded.
      */
     yInitialConditions.pop_front();
 };
diff --git a/ODElibrary/BDF.h b/ODElibrary/BDF.h
index 58ec50e..be2099d 100644
--- a/ODElibrary/BDF.h
+++ b/ODElibrary/BDF.h
@@ -1,33 +1,35 @@
 #ifndef ODELIBRARY_BDF_H
 #define ODELIBRARY_BDF_H
 
 #include <iostream>
 #include <Eigen/Dense>
 #include <list>
 #include <vector>
 #include "MultistepAbstractSolver.h"
 
 class BDF:public MultistepAbstractSolver{
 private:
     Eigen::MatrixXd updateMatrix;
 public:
     BDF();
     BDF(double t0,double tf,int n,int steps,Eigen::VectorXd& y0);
 
+    void setStepsAlgo(int n);
+    void setBCoeff();
     /**
      * This method overrides the setRightHandSide method of the abstract mother class. This operation is needed for implicit multistep
      * methods, because the matrix employed at each iteration for the update of y, updateMatrix, has to be computed and stored
      * as a variable of this class.
      * @param rhsMatrix: Matrix defining the ODE linear system
      * @param g_: function of t, part of the rhs expression.
      */
     void setRightHandSide(Eigen::MatrixXd & rhsMatrix, Eigen::VectorXd (*g_)(double t));
     /**
     * This function defines the update which is performed at each step of Adams Moulton method.
     * @param y : reference to the last computed value of y
     * @param t : time corresponding to the last computed value of y. The function computes y at time t+stepSize.
     */
     void step(Eigen::VectorXd &y, double t) override;
 };
 
 #endif
\ No newline at end of file
diff --git a/ODElibrary/ForwardEuler.cpp b/ODElibrary/ForwardEuler.cpp
index 8ab869e..ed32688 100644
--- a/ODElibrary/ForwardEuler.cpp
+++ b/ODElibrary/ForwardEuler.cpp
@@ -1,10 +1,17 @@
 
 #include <iostream>
 #include "ForwardEuler.h"
 
 
-
+ForwardEuler::ForwardEuler(double t0, double tf, int n, Eigen::VectorXd &y0){
+    initialTime = t0;
+    finalTime = tf;
+    nSteps = n;
+    stepSize = (tf - t0) / n;
+    dimension=y0.size();
+    initialValue = y0;
+}
 
 void ForwardEuler::step(Eigen::VectorXd &y, double t) {
     y += stepSize*(A*y + g(t)); // the simple Euler update
 }
diff --git a/ODElibrary/ForwardEuler.h b/ODElibrary/ForwardEuler.h
index e4ad8cb..c4ed5f1 100644
--- a/ODElibrary/ForwardEuler.h
+++ b/ODElibrary/ForwardEuler.h
@@ -1,20 +1,22 @@
 
 #ifndef ODELIBRARY_FORWARDEULER_H
 #define ODELIBRARY_FORWARDEULER_H
 #include "AbstractSolver.h"
 
 /**
  * Solver using the simple and very famous forward Euler method, defined as y_{t+1} = y_t +  h * (A*y_t + g(t))
  */
 class ForwardEuler : public AbstractSolver{
 public:
+    ForwardEuler():AbstractSolver(){};
+    ForwardEuler(double t0,double tf,int n,Eigen::VectorXd& y0);
     /**
      * Performs a step of the simple Euler method in place.
      * @param y : Eigen vector, the current value of the equation
      * @param t : double, the current time where the step must be performed
      */
     void step(Eigen::VectorXd &y, double t) override ;
 };
 
 #endif //ODELIBRARY_FORWARDEULER_H
 
diff --git a/ODElibrary/LinearCoefficients.cpp b/ODElibrary/LinearCoefficients.cpp
deleted file mode 100644
index bdd6650..0000000
--- a/ODElibrary/LinearCoefficients.cpp
+++ /dev/null
@@ -1,5 +0,0 @@
-//
-// Created by lionel on 25.11.19.
-//
-
-#include "LinearCoefficients.h"
diff --git a/ODElibrary/LinearCoefficients.h b/ODElibrary/LinearCoefficients.h
deleted file mode 100644
index 11a807f..0000000
--- a/ODElibrary/LinearCoefficients.h
+++ /dev/null
@@ -1,20 +0,0 @@
-//
-// Created by lionel on 25.11.19.
-//
-
-#ifndef ODELIBRARY_LINEARCOEFFICIENTS_H
-#define ODELIBRARY_LINEARCOEFFICIENTS_H
-
-
-class LinearCoefficients {
-private:
-    double ** A;
-    int dim;
-public:
-    LinearCoefficients();
-    ~LinearCoefficients();
-
-};
-
-
-#endif //ODELIBRARY_LINEARCOEFFICIENTS_H
diff --git a/ODElibrary/MultistepAbstractSolver.cpp b/ODElibrary/MultistepAbstractSolver.cpp
index f9a0609..8cb6670 100644
--- a/ODElibrary/MultistepAbstractSolver.cpp
+++ b/ODElibrary/MultistepAbstractSolver.cpp
@@ -1,148 +1,145 @@
 #include "MultistepAbstractSolver.h"
 
 MultistepAbstractSolver::MultistepAbstractSolver(){
     initialTime = 0;
     finalTime = 0;
     nSteps = 0;
     stepSize = 0;
     stepsAlgo=0;
     dimension=0;
+    convergenceRate=0;
 };
 
-MultistepAbstractSolver::MultistepAbstractSolver(double t0,double tf,int n,int steps,Eigen::VectorXd& y0) {
-    initialTime = t0;
-    finalTime = tf;
-    nSteps = n;
-    stepSize = (tf - t0) / n;
+MultistepAbstractSolver::MultistepAbstractSolver(double t0,double tf,int n,int steps,Eigen::VectorXd& y0):AbstractSolver(t0,tf,n,y0) {
     stepsAlgo=steps;
-    dimension=y0.size();
-    initialValue = y0;
     yInitialConditions.push_back(y0);
 };
 
-void MultistepAbstractSolver::setStepsAlgo(int n){
-    stepsAlgo=n;
-};
 
 int MultistepAbstractSolver::getStepsAlgo(){
     return stepsAlgo;
 };
 
 void MultistepAbstractSolver::computeMissingInitialValues(int orderRK){
     /**
      * This vector stores the coefficients needed by RK to compute the next value of y, yNext.
      * Such value is computed as yNext=yPrev+stepSize*(aCoeff[0]*k1+aCoeff[1]*k2+...), where
      * the ks correspond to evaluations of the rhs at some points of the interval (tPrev,tNext).
      * These points are some additional parameters of the Runge-Kutta method, chosen in order to guarantee
      * the desired convergence rate.
      * In order to make the remaining part of the code more compact, the vector aCoeff is always created with length 3,
      * as if needed for the highest order R-K method allowed by the code. Lower order R-K methods require fewer entries.
      * As a consequence, some of the entries are set equal to zero for lower order methods. This procedure implies
      * a few unneeded computations at line 95, but allows a more compact and easily extendable code.
      */
     Eigen::Vector4d aCoeff;
     switch(orderRK){
         case 0:
             {
+            rhsPreviousSteps.clear();
             Eigen::VectorXd rhsStart=A*initialValue+g(initialTime);
             rhsPreviousSteps.push_back(rhsStart);
             return;
             }
             break;
         case 1:
             aCoeff[0]=1, aCoeff[1]=0, aCoeff[2]=0,aCoeff[3]=0;
             break;
         case 2:
             aCoeff[0]=0.5, aCoeff[1]=0.5, aCoeff[2]=0,aCoeff[3]=0;
             break;
         case 3:
             aCoeff[0]=1./6., aCoeff[1]=4./6., aCoeff[2]=1./6.,aCoeff[3]=0;
             break;
         case 4:
             aCoeff[0]=1./4., aCoeff[1]=1./4., aCoeff[2]=1./4.,aCoeff[3]=1./4.;
             break;
         default:
             std::cerr<<"Order orderRK required too high. Try using a multistep method with fewer steps."<<"\n";
     }
     auto yPrev=initialValue;
     auto tPrev=initialTime;
     /**
      * \brief The evaluation of the rhs function at initialTime is computed and stored in rhsPreviousSteps.
+     * The list rhsPreviousSteps is cleared before starting to add vectors to it to be sure that no vector
+     * has already been added to it. This is necessary for the test performed in testAdamsBashforth.
      */
+    rhsPreviousSteps.clear();
     Eigen::VectorXd rhsStart=A*initialValue+g(initialTime);
     rhsPreviousSteps.push_back(rhsStart);
     /**
      * \brief Vectors k1, k2, k3, k4 define the updates of the R-K methods.
      */
     Eigen::VectorXd k1(dimension),k2(dimension),k3(dimension),k4(dimension);
     k2.setZero(), k3.setZero(),k4.setZero();
     int initialConditionsNeeded=orderRK;
     if(stepsAlgo==convergenceRate-1)
         initialConditionsNeeded-=1;
     for (int i=0;i<initialConditionsNeeded;i++){
         k1=A*yPrev+g(tPrev);
         /**
          * \brief In order to avoid computing the rhs function when unnecessary, this quantity is evaluated only if a higher order R-K method is actually needed.
          */
         if(orderRK>1)
             k2=A*(yPrev+0.5*stepSize*k1)+g(tPrev+0.5*stepSize); // With this choice, the order 2 method corresponds to Euler's modified method.
         if(orderRK>2)
             k3 = A * (yPrev + stepSize * (-1 * k1 + 2 * k2)) +g(tPrev + stepSize);// This choice for k3 allows us to "recycle" k1, k2, obtaining a 3rd order RK.
         if(orderRK>3) {
             k3=A *(yPrev+0.5*stepSize*k2)+g(tPrev+0.5*stepSize);
             k4 = A * (yPrev + stepSize*k3) + g(tPrev + stepSize);
         }
         Eigen::VectorXd yNext=yPrev+stepSize*(aCoeff[0]*k1+aCoeff[1]*k2+aCoeff[2]*k3+aCoeff[3]*k4);
         /**
          * \brief The value of the rhs function at tNext=tPrev+stepSize is stored in rhsPreviousSteps.
          */
         Eigen::VectorXd rhsNext=A*yNext+g(tPrev+stepSize);
         rhsPreviousSteps.push_back(rhsNext);
         /**
          * \brief The computed missing initial condition is stored in the list yInitialConditions.
          */
         yInitialConditions.push_back(yNext);
         yPrev=yNext;
         tPrev+=stepSize;
     }
 }
 
 void MultistepAbstractSolver::solve(std::ostream &stream){
+    yInitialConditions.clear();
+    yInitialConditions.push_back(initialValue);
     /**\brief
     * The missing initial values are computed, with a suitable one step algorithm.
     */
     if(stepsAlgo==convergenceRate-1)
         computeMissingInitialValues(convergenceRate-1);
     else
         computeMissingInitialValues(stepsAlgo-1);
     /**\brief
      * The initial conditions computed with computeMissingInitialValues are streamed to the output stream chosen by the user.
      */
     int count=0;
     auto it=yInitialConditions.begin();
     for(it=yInitialConditions.begin();it!=yInitialConditions.end();it++){
-        std::cout<<initialTime+stepSize*count<<" "<<(*it).transpose()<<"\n";
         stream<<initialTime+stepSize*count<<" "<<(*it).transpose()<<"\n";
         count++;
     }
     /**\brief
      * The variable tNext is the time at which y has to be computed at each step.
      * The variable yNext, instead, is initialized with the value of the last missing initial condition computed;
      * its value is progressively changed by the algorithm's iterations.
      */
     double t=initialTime+stepSize*(stepsAlgo-1);
     Eigen::VectorXd y = yInitialConditions.back();
     /** \brief
      * Since the first values of y have already been computed to obtain the needed initial conditions,
      * the integration algorithms has to be applied only (nSteps-stepsAlgo) times.
      */
     for (int i=0;i<=nSteps-stepsAlgo;i++) {
         step(y,t);
-        std::cout<<initialTime+stepSize*(i+stepsAlgo)<<" "<<y.transpose()<<"\n";
         stream<<initialTime+stepSize*(i+stepsAlgo)<<" "<<y.transpose()<<"\n";
+        t+=stepSize;
     }
 }
 
 void MultistepAbstractSolver::step(Eigen::VectorXd &y, double t) {
     throw std::runtime_error("function step not implemented in MultistepAbstractSolver: it's a purely virtual method, implemented in the daughter classes.");
 }
diff --git a/ODElibrary/MultistepAbstractSolver.h b/ODElibrary/MultistepAbstractSolver.h
index cd8e6a6..3e025cf 100644
--- a/ODElibrary/MultistepAbstractSolver.h
+++ b/ODElibrary/MultistepAbstractSolver.h
@@ -1,93 +1,90 @@
 /**
  * This abstract class is a daughter class of the AbstractSolver class.
  * It is  designed to be the mother class to its daughter classes AdamsBashforth, AdamsMoulton, BDF.
  * These classes will not be abstract and are the one that will be instanciated to solve ODE systems.
  * Creating an abstract multistep solver class allows to gather all the information and methods that are common
  * to all multistep algorithms.
  */
 
 
 
 #ifndef ODELIBRARY_MULTISTEPABSTRACTSOLVER_H
 #define ODELIBRARY_MULTISTEPABSTRACTSOLVER_H
 
 #include <iostream>
 #include <Eigen/Dense>
 #include <list>
 #include <vector>
 #include "AbstractSolver.h"
 
 
 class MultistepAbstractSolver:public AbstractSolver{
 protected:
     /**
      * This variable refers to the number of steps that characterizes the multistep algorithm.
      */
     int stepsAlgo;
     /**
      * This variable describes the convergence rate of the algorithm. It is of great importance since it provides
      * key information concerning the order of the method to be applied to compute the missing initial conditions;
      */
     int convergenceRate;
     /**
      * The vector bCoefficients stores the coefficients that determine the multistep algorithm. The coefficients
      * depend on the specific method chosen to solve the equation.
      */
     Eigen::VectorXd bCoefficients;
     /**
      * This list stores the y vectors that represent the initial conditions needed to perform
      * the first integration step of the algorithm.
      */
     std::list<Eigen::VectorXd> yInitialConditions;
     /**
      * This list stores only the values of the rhs function needed to perform the next integration step.
      * Its size depends on the number of steps of the multistep algorithm.
      */
      std::list<Eigen::VectorXd> rhsPreviousSteps;
 
 
 public:
     MultistepAbstractSolver();
     MultistepAbstractSolver(double t0,double tf,int n,int steps,Eigen::VectorXd& y0);
 
 
     /**\brief This function allows the user to check the number of steps of the multistep algorithm which is being used.
      *
      */
-
-    void setStepsAlgo(int n);
-
     int getStepsAlgo();
 
     /**
      * A multistep method of n steps needs n initial values to start its iterative process. Since only one initial value is usually
      * provided, i.e. the initial value corresponding to t=t0, the other missing values need to be computed. Such task
      * can be performed applying a one step method with convergence rate equal to the convergence rate of the multistep method minus 1.
      * By doing so, the numerical solution to the ODE system provided by the algorithm is guaranteed to converge to the true solution
      * with the convergence rate characterising the multistep method. In the context of this project, a suitable Runge-Kutta method is used.
      * @param orderRK: order of the RK method to be applied to compute missing initial values.
      */
     void computeMissingInitialValues(int orderRK);
 
     /**
      * This virtual method, which is then implemented by each daughter class in a specific way, solves the ODE system at hand.
      * @param stream: specifies where the computed values of the function have to be displayed/stored.
      */
-    void solve(std::ostream &stream = std::cout) override;
+    void solve(std::ostream &stream=std::cout) override;
     /**
      * This purely virtual method is implemented by the daughter classes AdamsBashforth, AdamsMoulton, BDF.
      * @param y : reference to the last computed value of y
      * @param t : time corresponding to the last computed value of y. The function computes y at time t+stepSize.
      */
     void step(Eigen::VectorXd &y, double t) override =0;
 
 };
 
 
 
 
 
 
 
 
 #endif
\ No newline at end of file
diff --git a/ODElibrary/RungeKutta.cpp b/ODElibrary/RungeKutta.cpp
index f80177e..c36aa7c 100644
--- a/ODElibrary/RungeKutta.cpp
+++ b/ODElibrary/RungeKutta.cpp
@@ -1,43 +1,53 @@
 
 #include <iostream>
 #include "RungeKutta.h"
 
+RungeKutta::RungeKutta(double t0,double tf,int n,int order, Eigen::VectorXd& y0){
+    initialTime = t0;
+    finalTime = tf;
+    nSteps = n;
+    stepSize = (tf - t0) / n;
+    dimension=y0.size();
+    initialValue = y0;
+    methodOrder=order;
+}
+
 void RungeKutta::step(Eigen::VectorXd &y, double t) {
     // The update is different depending on the odrer
     switch (methodOrder) {
         case 1 : // simple Euler step
             y += stepSize * (A * y + g(t));
             break;
         case 2 : // midpoint method
             y += stepSize * (A * (y + 0.5 * stepSize * (A * y + g(t))) + g(t + 0.5 * stepSize));
             break;
         case 3 : {// Kutta's third order method
             Eigen::VectorXd k1 = stepSize * (A * y + g(t));
             Eigen::VectorXd k2 = stepSize * (A * (y + k1 / 2.) + g(t + stepSize / 2.));
             Eigen::VectorXd k3 = stepSize * (A * (y - k1 + k2 / 2.) + g(t + stepSize));
             y += 1./6. * (k1 + 4 * k2 + k3);
             break;
         }
         case 4 : {//classical RK-4 method
             Eigen::VectorXd k1 = stepSize * (A * y + g(t));
             Eigen::VectorXd k2 = stepSize * (A * (y + k1 / 2.) + g(t + stepSize / 2.));
             Eigen::VectorXd k3 = stepSize * (A * (y + k2 / 2.) + g(t + stepSize / 2.));
             Eigen::VectorXd k4 = stepSize * (A * (y + k3) + g(t + stepSize));
             y += 1./6. * (k1 + 2 * k2 + 2 * k3 + k4);
             break;
         }
         default:
             throw std::domain_error("Runge-Kutta method order is not valid in step.");
     }
 }
 
 void RungeKutta::setMethodOrder(int n) {
     if (n < 1 || n > 4){ // check the bounds of implementation
         throw std::invalid_argument("Runge-Kutta method is only implemented for order from 1, 2, 3 or 4.");
     } else {
         methodOrder = n;
     }
 }
 int RungeKutta::getMethodOrder(){
     return methodOrder;
 };
diff --git a/ODElibrary/RungeKutta.h b/ODElibrary/RungeKutta.h
index 836d906..dd20c75 100644
--- a/ODElibrary/RungeKutta.h
+++ b/ODElibrary/RungeKutta.h
@@ -1,34 +1,35 @@
 
 #ifndef ODELIBRARY_RUNGEKUTTA_H
 #define ODELIBRARY_RUNGEKUTTA_H
 #include "AbstractSolver.h"
 
 /**
  * The famous Runge Kutta solver. Solutions for 4 methods order are available. It will compute the 4th order udpate by
  * default, a.k.a. the RK4 method.
  */
 class RungeKutta : public AbstractSolver{
 private:
     //
     int methodOrder;
 public:
+    RungeKutta(double t0,double tf,int n,int order, Eigen::VectorXd& y0);
     /**
      * Default constructor sets the method order to 4.
      */
     RungeKutta() : AbstractSolver() { methodOrder=4; };
     /**
      * Setter for the method order, will throw an exception if the order is invalid.
      * @param n : method order, integer, can be 1, 2, 3 or 4
      */
     void setMethodOrder(int n);
     /// getter for the method order
     int getMethodOrder();
     /**
      * Step method for runge Kutta, it will depend on the method order.
      * @param y : Eigen vector, the current value of the equation
      * @param t : double, the current time where the step must be performed
      */
     void step(Eigen::VectorXd &y, double t) override ;
 };
 
 #endif //ODELIBRARY_RUNGEKUTTA_H
\ No newline at end of file
diff --git a/tests/testAdamBashforth.cpp b/tests/testAdamBashforth.cpp
index 4389a1e..4fba515 100644
--- a/tests/testAdamBashforth.cpp
+++ b/tests/testAdamBashforth.cpp
@@ -1,112 +1,136 @@
-//
-// Created by lionel on 01.12.19.
-//
 
 # include <gtest/gtest.h>
 # include <cmath>
+#include <fstream>
 # include <cassert>
 # include <exception>
 #include "../ODElibrary/AdamsBashforth.h"
 
 // Function g(t) for testing
 Eigen::VectorXd ForcedOscillations(double t){
     Eigen::VectorXd s(2);
     s[0] = 0;
-    s[1] = cos(2*M_PI*t); // derivative of velocity
+    s[1] = sin(2*M_PI*t); // derivative of velocity
     return s;
 };
 
 // Fixture class
 class AdamsBashforthSolverTest : public ::testing::Test {
 protected:
     void SetUp() override {
     // create a solver
         solver.setInitialTime(-3);
         Eigen::VectorXd y0(2); // set the spring at rest
         y0[0] = 0;
         y0[1] = 0;
-        double k = 0.5; // spring constant
+        double k = 0.1; // spring constant
         Eigen::MatrixXd m(2,2); // from spring equation
         m(0,0)=0;
         m(0,1)=1;
         m(1,0)=-k;
         m(1,1)=0.3;
 
-        //AdamsBashforth solver=AdamsBashforth(0,20,2000,2,y0);
         solver.setInitialTime(0.);
-        solver.setFinalTime(20.);
-        solver.setStepSize(0.01);
-        solver.setRightHandSide(m, ForcedOscillations);
+        solver.setFinalTime(1.);
         solver.setInitialValue(y0);
-        solver.setStepsAlgo(2);
-
+        solver.setRightHandSide(m, &ForcedOscillations);
+        solver.setNumberOfSteps(100);
     }
 
 
     void TearDown() override {
         // clear particles
      //  solver;
     }
     AdamsBashforth solver;
+    double yFinal0Exact = 0.182532458;
+    double yFinal1Exact = 0.046017418;
 };
 
 
 // tests that we can assign and read a positive initial time
 TEST_F(AdamsBashforthSolverTest, initialAndFinalTimes) {
     solver.setInitialTime(0.5);
     EXPECT_EQ(solver.getInitialTime(), 0.5);
     solver.setFinalTime(4.2);
     EXPECT_EQ(solver.getFinalTime(), 4.2);
 }
 
 // tests that we can assign and read a positive timestep
 TEST_F(AdamsBashforthSolverTest, positiveTimestep) {
     solver.setStepSize(0.3);
     EXPECT_EQ(solver.getStepSize(), 0.3);
 }
 // tests we cannot assign a 0 or negative timestep
 TEST_F(AdamsBashforthSolverTest, nonPositiveTimestep) {
     ASSERT_THROW(solver.setStepSize(-0.1), std::domain_error);
     ASSERT_THROW(solver.setStepSize(0), std::domain_error);
 }
 
 // tests we cannot assign a different size inputs for the initial value
 TEST_F(AdamsBashforthSolverTest, differentSizedInitialValue) {
     Eigen::VectorXd yTooLong(3);
     ASSERT_THROW(solver.setInitialValue(yTooLong), std::length_error);
 }
 
 // tests we cannot assign a different size inputs for the rhs matrix
 TEST_F(AdamsBashforthSolverTest, differentSizedRHS) {
     Eigen::MatrixXd mFalse1(2,3);
     Eigen::MatrixXd mFalse2(1,2);
     Eigen::MatrixXd mFalse3(3,3);
     ASSERT_THROW(solver.setRightHandSide(mFalse1, ForcedOscillations), std::length_error);
     ASSERT_THROW(solver.setRightHandSide(mFalse2, ForcedOscillations), std::length_error);
     ASSERT_THROW(solver.setRightHandSide(mFalse3, ForcedOscillations), std::length_error);
 }
 
 // tests we cannot assign a rhs function
 Eigen::VectorXd RHSgTooLong(double t){
     Eigen::VectorXd s(4);
     return s;
 };
 Eigen::VectorXd RHSgThrowingException(double t){
     Eigen::VectorXd s(2);
     throw std::exception();
     return s;
 };
 Eigen::VectorXd RHSgThrowingInteger(double t){
     Eigen::VectorXd s(2);
     throw 42;
     return s;
 };
 TEST_F(AdamsBashforthSolverTest, timeDependentRHSFunction) {
     Eigen::MatrixXd m(2,2);
 
     ASSERT_THROW(solver.setRightHandSide(m, RHSgTooLong), std::length_error);
     ASSERT_THROW(solver.setRightHandSide(m, RHSgThrowingException), std::runtime_error);
     ASSERT_THROW(solver.setRightHandSide(m, RHSgThrowingInteger), std::runtime_error);
 
 }
 
+
+TEST_F(AdamsBashforthSolverTest, testSolve) {
+    for(int s=1;s<=4;s++) {
+        solver.setStepsAlgo(s);
+        // writes the output to the file
+        std::ofstream testFile;
+        testFile.open("ABSolverTest.dat");
+        assert(testFile.is_open());
+        solver.solve(testFile);
+        testFile.close();
+        // read the file where the output is
+        std::ifstream readFile("ABSolverTest.dat");
+        assert(readFile.is_open());
+        double t = 0, x = 0, y = 0;
+        int i = 0;
+        while (!readFile.eof()) {
+            readFile >> t >> x >> y;
+            i++;
+        }
+        readFile.close();
+        // check the value are close enough to the specified values
+        //std::cout << t << " " << x << " " << y << " " << i;
+        ASSERT_NEAR(x, yFinal0Exact, 1e-2);
+        ASSERT_NEAR(y, yFinal1Exact, 1e-2);
+    }
+
+}
diff --git a/tests/testAdamsMoulton.cpp b/tests/testAdamsMoulton.cpp
index d65989a..416abbf 100644
--- a/tests/testAdamsMoulton.cpp
+++ b/tests/testAdamsMoulton.cpp
@@ -1,106 +1,143 @@
 # include <gtest/gtest.h>
 # include <cmath>
 # include <cassert>
+#include <fstream>
 # include <exception>
 #include "../ODElibrary/AdamsMoulton.h"
 
 // Function g(t) for testing
 Eigen::VectorXd ForcedOscillations(double t){
     Eigen::VectorXd s(2);
     s[0] = 0;
-    s[1] = cos(2*M_PI*t); // derivative of velocity
+    s[1] = sin(2*M_PI*t); // derivative of velocity
     return s;
 };
 
 // Fixture class
 class AdamsMoultonSolverTest : public ::testing::Test {
 protected:
     void SetUp() override {
         // create a solver
         solver.setInitialTime(-3);
         Eigen::VectorXd y0(2); // set the spring at rest
         y0[0] = 0;
         y0[1] = 0;
-        double k = 0.5; // spring constant
+        double k = 0.1; // spring constant
         Eigen::MatrixXd m(2,2); // from spring equation
         m(0,0)=0;
         m(0,1)=1;
         m(1,0)=-k;
         m(1,1)=0.3;
 
         solver.setInitialTime(0.);
-        solver.setFinalTime(20.);
-        solver.setStepSize(0.01);
+        solver.setFinalTime(1.);
+        solver.setNumberOfSteps(100);
+        solver.setStepsAlgo(2);
         solver.setRightHandSide(m, ForcedOscillations);
         solver.setInitialValue(y0);
-        solver.setStepsAlgo(2);
 
     }
 
 
     void TearDown() override {
         // clear particles
         //  solver;
     }
     AdamsMoulton solver;
+    double yFinal0Exact = 0.182532458;
+    double yFinal1Exact = 0.046017418;
 };
 
 
 // tests that we can assign and read a positive initial time
 TEST_F(AdamsMoultonSolverTest, initialAndFinalTimes) {
     solver.setInitialTime(0.5);
     EXPECT_EQ(solver.getInitialTime(), 0.5);
     solver.setFinalTime(4.2);
     EXPECT_EQ(solver.getFinalTime(), 4.2);
 }
 
 // tests that we can assign and read a positive timestep
 TEST_F(AdamsMoultonSolverTest, positiveTimestep) {
     solver.setStepSize(0.3);
     EXPECT_EQ(solver.getStepSize(), 0.3);
 }
 // tests we cannot assign a 0 or negative timestep
 TEST_F(AdamsMoultonSolverTest, nonPositiveTimestep) {
     ASSERT_THROW(solver.setStepSize(-0.1), std::domain_error);
     ASSERT_THROW(solver.setStepSize(0), std::domain_error);
 }
 
 // tests we cannot assign a different size inputs for the initial value
 TEST_F(AdamsMoultonSolverTest, differentSizedInitialValue) {
     Eigen::VectorXd yTooLong(3);
     ASSERT_THROW(solver.setInitialValue(yTooLong), std::length_error);
 }
 
 // tests we cannot assign a different size inputs for the rhs matrix
 TEST_F(AdamsMoultonSolverTest, differentSizedRHS) {
     Eigen::MatrixXd mFalse1(2,3);
     Eigen::MatrixXd mFalse2(1,2);
     Eigen::MatrixXd mFalse3(3,3);
     ASSERT_THROW(solver.setRightHandSide(mFalse1, ForcedOscillations), std::length_error);
     ASSERT_THROW(solver.setRightHandSide(mFalse2, ForcedOscillations), std::length_error);
     ASSERT_THROW(solver.setRightHandSide(mFalse3, ForcedOscillations), std::length_error);
 }
 
 // tests we cannot assign a rhs function
 Eigen::VectorXd RHSgTooLong(double t){
     Eigen::VectorXd s(4);
     return s;
 };
 Eigen::VectorXd RHSgThrowingException(double t){
     Eigen::VectorXd s(2);
     throw std::exception();
     return s;
 };
 Eigen::VectorXd RHSgThrowingInteger(double t){
     Eigen::VectorXd s(2);
     throw 42;
     return s;
 };
 TEST_F(AdamsMoultonSolverTest, timeDependentRHSFunction) {
     Eigen::MatrixXd m(2,2);
 
     ASSERT_THROW(solver.setRightHandSide(m, RHSgTooLong), std::length_error);
     ASSERT_THROW(solver.setRightHandSide(m, RHSgThrowingException), std::runtime_error);
     ASSERT_THROW(solver.setRightHandSide(m, RHSgThrowingInteger), std::runtime_error);
 
 }
+
+
+TEST_F(AdamsMoultonSolverTest, testSolve) {
+    for(int s=1;s<=4;s++) {
+        solver.setStepsAlgo(s);
+        Eigen::MatrixXd m(2,2); // from spring equation
+        m(0,0)=0;
+        m(0,1)=1;
+        m(1,0)=-0.1;
+        m(1,1)=0.3;
+        solver.setRightHandSide(m, ForcedOscillations);
+        // writes the output to the file
+        std::ofstream testFile;
+        testFile.open("AMSolverTest.dat");
+        assert(testFile.is_open());
+        solver.solve(testFile);
+        testFile.close();
+        // read the file where the output is
+        std::ifstream readFile("AMSolverTest.dat");
+        assert(readFile.is_open());
+        double t = 0, x = 0, y = 0;
+        int i = 0;
+        while (!readFile.eof()) {
+            readFile >> t >> x >> y;
+            i++;
+        }
+        readFile.close();
+        // check the value are close enough to the specified values
+        //std::cout << t << " " << x << " " << y << " " << i;
+        ASSERT_NEAR(x, yFinal0Exact, 2e-2);
+        ASSERT_NEAR(y, yFinal1Exact, 2e-2);
+    }
+
+}
diff --git a/tests/testBDF.cpp b/tests/testBDF.cpp
index 606b9d1..eec88db 100644
--- a/tests/testBDF.cpp
+++ b/tests/testBDF.cpp
@@ -1,106 +1,142 @@
 # include <gtest/gtest.h>
 # include <cmath>
 # include <cassert>
 # include <exception>
+#include <fstream>
 #include "../ODElibrary/BDF.h"
 
 // Function g(t) for testing
 Eigen::VectorXd ForcedOscillations(double t){
     Eigen::VectorXd s(2);
     s[0] = 0;
-    s[1] = cos(2*M_PI*t); // derivative of velocity
+    s[1] = sin(2*M_PI*t); // derivative of velocity
     return s;
 };
 
 // Fixture class
 class BDFSolverTest : public ::testing::Test {
 protected:
     void SetUp() override {
         // create a solver
         solver.setInitialTime(-3);
         Eigen::VectorXd y0(2); // set the spring at rest
         y0[0] = 0;
         y0[1] = 0;
-        double k = 0.5; // spring constant
+        double k = 0.1; // spring constant
         Eigen::MatrixXd m(2,2); // from spring equation
         m(0,0)=0;
         m(0,1)=1;
         m(1,0)=-k;
         m(1,1)=0.3;
 
         solver.setInitialTime(0.);
-        solver.setFinalTime(20.);
-        solver.setStepSize(0.01);
+        solver.setFinalTime(1.);
+        solver.setNumberOfSteps(100);
+        solver.setStepsAlgo(2);
         solver.setRightHandSide(m, ForcedOscillations);
         solver.setInitialValue(y0);
-        solver.setStepsAlgo(2);
 
     }
 
 
     void TearDown() override {
         // clear particles
         //  solver;
     }
     BDF solver;
+    double yFinal0Exact = 0.182532458;
+    double yFinal1Exact = 0.046017418;
 };
 
 
 // tests that we can assign and read a positive initial time
 TEST_F(BDFSolverTest, initialAndFinalTimes) {
     solver.setInitialTime(0.5);
     EXPECT_EQ(solver.getInitialTime(), 0.5);
     solver.setFinalTime(4.2);
     EXPECT_EQ(solver.getFinalTime(), 4.2);
 }
 
 // tests that we can assign and read a positive timestep
 TEST_F(BDFSolverTest, positiveTimestep) {
     solver.setStepSize(0.3);
     EXPECT_EQ(solver.getStepSize(), 0.3);
 }
 // tests we cannot assign a 0 or negative timestep
 TEST_F(BDFSolverTest, nonPositiveTimestep) {
     ASSERT_THROW(solver.setStepSize(-0.1), std::domain_error);
     ASSERT_THROW(solver.setStepSize(0), std::domain_error);
 }
 
 // tests we cannot assign a different size inputs for the initial value
 TEST_F(BDFSolverTest, differentSizedInitialValue) {
     Eigen::VectorXd yTooLong(3);
     ASSERT_THROW(solver.setInitialValue(yTooLong), std::length_error);
 }
 
 // tests we cannot assign a different size inputs for the rhs matrix
 TEST_F(BDFSolverTest, differentSizedRHS) {
     Eigen::MatrixXd mFalse1(2,3);
     Eigen::MatrixXd mFalse2(1,2);
     Eigen::MatrixXd mFalse3(3,3);
     ASSERT_THROW(solver.setRightHandSide(mFalse1, ForcedOscillations), std::length_error);
     ASSERT_THROW(solver.setRightHandSide(mFalse2, ForcedOscillations), std::length_error);
     ASSERT_THROW(solver.setRightHandSide(mFalse3, ForcedOscillations), std::length_error);
 }
 
 // tests we cannot assign a rhs function
 Eigen::VectorXd RHSgTooLong(double t){
     Eigen::VectorXd s(4);
     return s;
 };
 Eigen::VectorXd RHSgThrowingException(double t){
     Eigen::VectorXd s(2);
     throw std::exception();
     return s;
 };
 Eigen::VectorXd RHSgThrowingInteger(double t){
     Eigen::VectorXd s(2);
     throw 42;
     return s;
 };
 TEST_F(BDFSolverTest, timeDependentRHSFunction) {
     Eigen::MatrixXd m(2,2);
 
     ASSERT_THROW(solver.setRightHandSide(m, RHSgTooLong), std::length_error);
     ASSERT_THROW(solver.setRightHandSide(m, RHSgThrowingException), std::runtime_error);
     ASSERT_THROW(solver.setRightHandSide(m, RHSgThrowingInteger), std::runtime_error);
 
-}
\ No newline at end of file
+}
+
+TEST_F(BDFSolverTest, testSolve) {
+    for(int s=1;s<=5;s++) {
+        solver.setStepsAlgo(s);
+        Eigen::MatrixXd m(2,2); // from spring equation
+        m(0,0)=0;
+        m(0,1)=1;
+        m(1,0)=-0.1;
+        m(1,1)=0.3;
+        solver.setRightHandSide(m, ForcedOscillations);
+        // writes the output to the file
+        std::ofstream testFile;
+        testFile.open("BDFSolverTest.dat");
+        assert(testFile.is_open());
+        solver.solve(testFile);
+        testFile.close();
+        // read the file where the output is
+        std::ifstream readFile("BDFSolverTest.dat");
+        assert(readFile.is_open());
+        double t = 0, x = 0, y = 0;
+        int i = 0;
+        while (!readFile.eof()) {
+            readFile >> t >> x >> y;
+            i++;
+        }
+        readFile.close();
+        // check the value are close enough to the specified values
+        //std::cout << t << " " << x << " " << y << " " << i;
+        ASSERT_NEAR(x, yFinal0Exact, 2e-2);
+        ASSERT_NEAR(y, yFinal1Exact, 2e-2);
+    }
+
+}