diff --git a/ODElibrary/AdamsBashforth.cpp b/ODElibrary/AdamsBashforth.cpp
index 935ec4c..c651ba6 100644
--- a/ODElibrary/AdamsBashforth.cpp
+++ b/ODElibrary/AdamsBashforth.cpp
@@ -1,48 +1,86 @@
 #include "AdamsBashforth.h"
 //#include "MultistepAbstractSolver.h"
 
 AdamsBashforth::AdamsBashforth(double t0,double tf,int n,int steps,Eigen::VectorXd& y0):MultistepAbstractSolver(t0,tf,n,steps,y0){
     /**\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){
         ///When steps=1, Adams Bashforth multistep algorithm actually coincides with Forward Euler.
         case 1:
             bCoefficients[0]=1;
             break;
         case 2:
             bCoefficients[0]=-0.5,bCoefficients[1]=1.5;
             break;
         case 3:
-            bCoefficients[0]=25./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.;
             break;
         default:
             std::cerr<<"The number of steps has to be an integer between 1 and 4."<<"\n";
     }
-    /**\brief
-     * The missing intial values are computed, with a suitable one step algorithm.
-     */
-    computeMissingInitialValues(steps-1);
-}
+};
 
 void AdamsBashforth::solve(std::ostream &stream){
-    double tNext=initialTime+stepSize*stepsAlgo;
-    auto yNext = *(*yPreviousSteps.end());
+    /**\brief
+     * The missing initial values are computed, with a suitable one step algorithm.
+     */
+    computeMissingInitialValues(stepsAlgo-1);
+    /**\brief
+     * The initial conditions computed with computeMissingInitialValues are streamed to the output stream chosen by the user.
+     */
     int count=0;
-    for (int step=0;step<nSteps-stepsAlgo;step++) {
-        auto it=rhsPreviousSteps.begin();
+    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 tNext=initialTime+stepSize*stepsAlgo;
+    Eigen::VectorXd yNext = *(yInitialConditions.end()-1);
+    /** \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 step=0;step<=nSteps-stepsAlgo;step++) {
+        auto it2=rhsPreviousSteps.begin();
         count=0;
-        for (it=rhsPreviousSteps.begin();it!=rhsPreviousSteps.end();it++) {
-            yNext +=(*(*it))*stepSize*bCoefficients[count]; //reversed coefficients in bCoefficients!
+        /**\brief
+         * The following cycle, which coincides with one integration step, defines the update of y.
+         */
+        for (it2=rhsPreviousSteps.begin();it2!=rhsPreviousSteps.end();it2++) {
+            yNext +=(*it2)*stepSize*bCoefficients[count]; //reversed coefficients in bCoefficients!
             count++;
         }
+        /**\brief
+         * The computed value is then streamed to the output stream selected by the user.
+         */
+        std::cout<<initialTime+stepSize*(step+stepsAlgo)<<" "<<yNext.transpose()<<"\n";
+        stream<<initialTime+stepSize*(step+stepsAlgo)<<" "<<yNext.transpose()<<"\n";
+        /**\brief
+         * Having computed yNext, the rhs function can subsequently be evaluated at tNext.
+         */
         Eigen::VectorXd rhsNext = A *yNext+g(tNext);//+g();
-        rhsPreviousSteps.push_back(&rhsNext);
-        rhsPreviousSteps.erase(rhsPreviousSteps.begin());
+        /**\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();
         tNext+=stepSize;
     }
 }
+
+void AdamsBashforth::solve(){
+
+};
diff --git a/ODElibrary/AdamsBashforth.h b/ODElibrary/AdamsBashforth.h
index 56f8aa1..dd3f8af 100644
--- a/ODElibrary/AdamsBashforth.h
+++ b/ODElibrary/AdamsBashforth.h
@@ -1,20 +1,21 @@
 #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(double t0,double tf,int n,int steps,Eigen::VectorXd& y0);
 
+    void solve() override;
     void solve(std::ostream &stream) override;
 };
 
 
 
 #endif
\ No newline at end of file
diff --git a/ODElibrary/MultistepAbstractSolver.cpp b/ODElibrary/MultistepAbstractSolver.cpp
index 1209f94..3faf7b8 100644
--- a/ODElibrary/MultistepAbstractSolver.cpp
+++ b/ODElibrary/MultistepAbstractSolver.cpp
@@ -1,68 +1,79 @@
 #include "MultistepAbstractSolver.h"
 
 MultistepAbstractSolver::MultistepAbstractSolver(double t0,double tf,int n,int steps,Eigen::VectorXd& y0) {
     initialTime = t0;
     finalTime = tf;
     nSteps = n;
     stepSize = (tf - t0) / n;
     stepsAlgo=steps;
-    Eigen::VectorXd v(steps);
-    bCoefficients = v;
+    dimension=y0.size();
     initialValue = y0;
-    yPreviousSteps.push_back(&y0); ///is it right? ask.
+    yInitialConditions.push_back(y0); ///is it right? ask.
+};
+
+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 order.
+     * 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 may require fewer entries.
-     * As a consequence, some of the entries are equal to zero for lower order methods. This procedure implies
-     * a few unneeded computations but allows a more compact and easily extendable.
+     * 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 68, but allows a more compact and easily extendable code.
      */
     Eigen::Vector3d aCoeff;
     switch(orderRK){
         case 1:
             aCoeff[0]=1, aCoeff[1]=0, aCoeff[2]=0;
             break;
         case 2:
             aCoeff[0]=0.5, aCoeff[1]=0.5, aCoeff[2]=0;
             break;
         case 3:
             aCoeff[0]=(double)1/6, aCoeff[1]=(double)4/6, aCoeff[2]=(double)1/6;
             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 A pointer to the value of the rhs function at initialTime is stored in rhsPreviousSteps.
+     * \brief The evaluation of the rhs function at initialTime is computed and stored in rhsPreviousSteps.
      */
     Eigen::VectorXd rhsStart=A*initialValue+g(initialTime);
-    rhsPreviousSteps.push_back(&rhsStart);
+    rhsPreviousSteps.push_back(rhsStart);
     /**
-     * \brief Vectors k2, k3 define the updates of the higher order R-K methods.
+     * \brief Vectors k1, k2, k3 define the updates of the R-K methods.
      */
-    Eigen::VectorXd k2(dimension),k3(dimension);
+    Eigen::VectorXd k1(dimension),k2(dimension),k3(dimension);
     k2.setZero(), k3.setZero();
     for (int i=0;i<orderRK;i++){
-        auto k1=A*yPrev+g(tPrev);
+        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)
-            auto k2=A*(yPrev+0.5*stepSize*k1)+g(tPrev+0.5*stepSize); // With this choice, the order two method corresponds to Euler's modified method.
+            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)
-            auto k3=A*(yPrev+stepSize*(-1*k1+2*k2))+g(tPrev+stepSize);// This choice for k3 allows us to "recycle" k1, k2, obtaining a 3rd oredr RK.
+            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.
+
         Eigen::VectorXd yNext=yPrev+stepSize*(aCoeff[0]*k1+aCoeff[1]*k2+aCoeff[2]*k3);  // why 'auto' does not work? ->check if it is right
         /**
          * \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);
-        //yPreviousSteps.push_back(&yNext);
+        rhsPreviousSteps.push_back(rhsNext);
+        /**
+         * \brief The computed missing initial condition is stored in the list yInitialConditions.
+         */
+        yInitialConditions.push_back(yNext);
+        yPrev=yNext;
         tPrev+=stepSize;
     }
 }
\ No newline at end of file
diff --git a/ODElibrary/MultistepAbstractSolver.h b/ODElibrary/MultistepAbstractSolver.h
index 7d1c477..4ef85c7 100644
--- a/ODElibrary/MultistepAbstractSolver.h
+++ b/ODElibrary/MultistepAbstractSolver.h
@@ -1,71 +1,80 @@
-
 /**
- * This abstract class is daughter of the AbstractSolver class.
+ * 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;
     /**
-     * bCoefficients stores the coefficients that determine the multistep algorithm. The coefficients depend on the
-     * specific method chosen to solve the equation.
+     * 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 vector stores pointers to the values of y needed to perform the next integration step.
-     * The size of this vector depends on the number of steps that characterizes the multistep algorithm.
-     */
-    std::vector<Eigen::VectorXd*> yPreviousSteps;
-    /**
-     * This vector stores pointers to the values of the rhs function needed to perform the next integration step.
-     * Its size depends on the number of steps of the algorithm.
+     * This vector stores the y vectors that represent the initial conditions needed to perform
+     * the first integration step of the algorithm. A vector is used instead of a list to make
+     * it easier to access the last element of it, which is a needed operation in the algorithm.
+     * (A reverse iterator on a list could have also been used.)
      */
-     std::vector<Eigen::VectorXd*> rhsPreviousSteps;
+    std::vector<Eigen::VectorXd> yInitialConditions;  //try
     /**
-     * This vector stores the values of the vectorial function g(t) needed to perform the next integration step.
-     * The size depends on the number of steps that characterizes the multistep algorithm.
+     * 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::vector<Eigen::VectorXd*> gPreviousSteps;
+     std::list<Eigen::VectorXd> rhsPreviousSteps;
+
 
 public:
     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.
+     *
+     */
+    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.
+     */
     virtual void solve(std::ostream &stream) = 0;
 
 };
 
 
 
 
 
 
 
 
 #endif
\ No newline at end of file