diff --git a/ModelEvaluatorImpl.hpp b/ModelEvaluatorImpl.hpp
index 3e9a4cef6..0f8969f0b 100644
--- a/ModelEvaluatorImpl.hpp
+++ b/ModelEvaluatorImpl.hpp
@@ -77,28 +77,27 @@ namespace ModelLib
         typedef typename ModelEvaluator<ScalarT, IdxT>::real_type real_type;
 
         ModelEvaluatorImpl()
-            : size_(0),
-              size_quad_(0),
-              size_opt_(0)
-        {
-        }
+          : size_(0),
+            size_quad_(0),
+            size_opt_(0)
+        {}
 
         ModelEvaluatorImpl(IdxT size, IdxT size_quad, IdxT size_opt)
-            : size_(size),
-              size_quad_(size_quad),
-              size_opt_(size_opt),
-              y_(size_),
-              yp_(size_),
-              f_(size_),
-              g_(size_quad_),
-              yB_(size_),
-              ypB_(size_),
-              fB_(size_),
-              gB_(size_opt_),
-              J_(COO_Matrix<ScalarT, IdxT>()),
-              param_(size_opt_),
-              param_up_(size_opt_),
-              param_lo_(size_opt_)
+          : size_(size),
+            size_quad_(size_quad),
+            size_opt_(size_opt),
+            y_(size_),
+            yp_(size_),
+            f_(size_),
+            g_(size_quad_),
+            yB_(size_),
+            ypB_(size_),
+            fB_(size_),
+            gB_(size_opt_),
+            J_(COO_Matrix<ScalarT, IdxT>()),
+            param_(size_opt_),
+            param_up_(size_opt_),
+            param_lo_(size_opt_)
         {
         }
 
@@ -134,143 +133,143 @@ namespace ModelLib
         //     std::cout << "updateTime: t = " << time_ << "\n";
         // }
 
-        virtual void setTolerances(real_type &rtol, real_type &atol) const
+        virtual void setTolerances(real_type& rtol, real_type& atol) const
         {
             rtol = rtol_;
             atol = atol_;
         }
 
-        virtual void setMaxSteps(IdxT &msa) const
+        virtual void setMaxSteps(IdxT& msa) const
         {
             msa = max_steps_;
         }
 
-        std::vector<ScalarT> &y()
+        std::vector<ScalarT>& y()
         {
             return y_;
         }
 
-        const std::vector<ScalarT> &y() const
+        const std::vector<ScalarT>& y() const
         {
             return y_;
         }
 
-        std::vector<ScalarT> &yp()
+        std::vector<ScalarT>& yp()
         {
             return yp_;
         }
 
-        const std::vector<ScalarT> &yp() const
+        const std::vector<ScalarT>& yp() const
         {
             return yp_;
         }
 
-        std::vector<bool> &tag()
+        std::vector<bool>& tag()
         {
             return tag_;
         }
 
-        const std::vector<bool> &tag() const
+        const std::vector<bool>& tag() const
         {
             return tag_;
         }
 
-        std::vector<ScalarT> &yB()
+        std::vector<ScalarT>& yB()
         {
             return yB_;
         }
 
-        const std::vector<ScalarT> &yB() const
+        const std::vector<ScalarT>& yB() const
         {
             return yB_;
         }
 
-        std::vector<ScalarT> &ypB()
+        std::vector<ScalarT>& ypB()
         {
             return ypB_;
         }
 
-        const std::vector<ScalarT> &ypB() const
+        const std::vector<ScalarT>& ypB() const
         {
             return ypB_;
         }
 
-        std::vector<ScalarT> &param()
+        std::vector<ScalarT>& param()
         {
             return param_;
         }
 
-        const std::vector<ScalarT> &param() const
+        const std::vector<ScalarT>& param() const
         {
             return param_;
         }
 
-        std::vector<ScalarT> &param_up()
+        std::vector<ScalarT>& param_up()
         {
             return param_up_;
         }
 
-        const std::vector<ScalarT> &param_up() const
+        const std::vector<ScalarT>& param_up() const
         {
             return param_up_;
         }
 
-        std::vector<ScalarT> &param_lo()
+        std::vector<ScalarT>& param_lo()
         {
             return param_lo_;
         }
 
-        const std::vector<ScalarT> &param_lo() const
+        const std::vector<ScalarT>& param_lo() const
         {
             return param_lo_;
         }
 
-        std::vector<ScalarT> &getResidual()
+        std::vector<ScalarT>& getResidual()
         {
             return f_;
         }
 
-        const std::vector<ScalarT> &getResidual() const
+        const std::vector<ScalarT>& getResidual() const
         {
             return f_;
         }
 
-        COO_Matrix<ScalarT, IdxT> &getJacobian()
+        COO_Matrix<ScalarT, IdxT>& getJacobian()
         {
             return J_;
         }
 
-        const COO_Matrix<ScalarT, IdxT> &getJacobian() const
+        const COO_Matrix<ScalarT, IdxT>& getJacobian() const
         {
             return J_;
         }
 
-        std::vector<ScalarT> &getIntegrand()
+        std::vector<ScalarT>& getIntegrand()
         {
             return g_;
         }
 
-        const std::vector<ScalarT> &getIntegrand() const
+        const std::vector<ScalarT>& getIntegrand() const
         {
             return g_;
         }
 
-        std::vector<ScalarT> &getAdjointResidual()
+        std::vector<ScalarT>& getAdjointResidual()
         {
             return fB_;
         }
 
-        const std::vector<ScalarT> &getAdjointResidual() const
+        const std::vector<ScalarT>& getAdjointResidual() const
         {
             return fB_;
         }
 
-        std::vector<ScalarT> &getAdjointIntegrand()
+        std::vector<ScalarT>& getAdjointIntegrand()
         {
             return gB_;
         }
 
-        const std::vector<ScalarT> &getAdjointIntegrand() const
+        const std::vector<ScalarT>& getAdjointIntegrand() const
         {
             return gB_;
         }
@@ -281,6 +280,8 @@ namespace ModelLib
             return idc_;
         }
 
+
+
     protected:
         IdxT size_;
         IdxT nnz_;
@@ -313,8 +314,10 @@ namespace ModelLib
         IdxT max_steps_;
 
         IdxT idc_;
+
     };
 
+
 } // namespace ModelLib
 
 #endif // _MODEL_EVALUATOR_IMPL_HPP_
\ No newline at end of file
diff --git a/Solver/Dynamic/Ida.cpp b/Solver/Dynamic/Ida.cpp
index 3ae6854c4..fbe7c8973 100644
--- a/Solver/Dynamic/Ida.cpp
+++ b/Solver/Dynamic/Ida.cpp
@@ -57,703 +57,718 @@
  *
  */
 
+
 #include <iostream>
 #include <iomanip>
 
-#include <idas/idas_direct.h> /* access to IDADls interface           */
+#include <idas/idas_direct.h>              /* access to IDADls interface           */
 #include <idas/idas.h>
 
-// Sundials Sparse KLU
+//Sundials Sparse KLU
 #include <sunmatrix/sunmatrix_sparse.h>
-#include <sunlinsol/sunlinsol_klu.h>
+#include <sunlinsol/sunlinsol_klu.h>   
 
 #include "ModelEvaluator.hpp"
 #include "Ida.hpp"
 
+
 namespace AnalysisManager
 {
 
-    namespace Sundials
-    {
+namespace Sundials
+{
 
-        template <class ScalarT, typename IdxT>
-        Ida<ScalarT, IdxT>::Ida(ModelLib::ModelEvaluator<ScalarT, IdxT> *model) : DynamicSolver<ScalarT, IdxT>(model)
-        {
-            int retval = 0;
+    template <class ScalarT, typename IdxT>
+    Ida<ScalarT, IdxT>::Ida(ModelLib::ModelEvaluator<ScalarT, IdxT>* model) : DynamicSolver<ScalarT, IdxT>(model)
+    {
+        int retval = 0;
 
-            // Create the SUNDIALS context that all SUNDIALS objects require
-            retval = SUNContext_Create(NULL, &context_);
-            checkOutput(retval, "SUNContext");
-            solver_ = IDACreate(context_);
-            tag_ = NULL;
-        }
+        // Create the SUNDIALS context that all SUNDIALS objects require
+        retval = SUNContext_Create(NULL, &context_);
+        checkOutput(retval, "SUNContext");
+        solver_ = IDACreate(context_);
+        tag_ = NULL;
+    }
 
-        template <class ScalarT, typename IdxT>
-        Ida<ScalarT, IdxT>::~Ida()
-        {
-        }
+    template <class ScalarT, typename IdxT>
+    Ida<ScalarT, IdxT>::~Ida()
+    {
+    }
 
-        template <class ScalarT, typename IdxT>
-        int Ida<ScalarT, IdxT>::configureSimulation()
-        {
-            int retval = 0;
+    template <class ScalarT, typename IdxT>
+    int Ida<ScalarT, IdxT>::configureSimulation()
+    {
+        int retval = 0;
 
-            // Allocate solution vectors
-            yy_ = N_VNew_Serial(model_->size(), context_);
-            checkAllocation((void *)yy_, "N_VNew_Serial");
-            yp_ = N_VClone(yy_);
-            checkAllocation((void *)yp_, "N_VClone");
+        // Allocate solution vectors
+        yy_ = N_VNew_Serial(model_->size(), context_);
+        checkAllocation((void*) yy_, "N_VNew_Serial");
+        yp_ = N_VClone(yy_);
+        checkAllocation((void*) yp_, "N_VClone");
 
-            // get intial conditions
-            this->getDefaultInitialCondition();
+        //get intial conditions
+        this->getDefaultInitialCondition();
 
-            // Create vectors to store restart initial condition
-            yy0_ = N_VClone(yy_);
-            checkAllocation((void *)yy0_, "N_VClone");
-            yp0_ = N_VClone(yy_);
-            checkAllocation((void *)yp0_, "N_VClone");
+        // Create vectors to store restart initial condition
+        yy0_ = N_VClone(yy_);
+        checkAllocation((void*) yy0_, "N_VClone");
+        yp0_ = N_VClone(yy_);
+        checkAllocation((void*) yp0_, "N_VClone");
 
-            // Dummy initial time; will be overridden.
-            const realtype t0 = RCONST(0.0);
+        // Dummy initial time; will be overridden.
+        const realtype t0 = RCONST(0.0);
 
-            // Allocate and initialize IDA workspace
-            retval = IDAInit(solver_, this->Residual, t0, yy_, yp_);
-            checkOutput(retval, "IDAInit");
+        // Allocate and initialize IDA workspace
+        retval = IDAInit(solver_, this->Residual, t0, yy_, yp_);
+        checkOutput(retval, "IDAInit");
 
-            // Set pointer to model data
-            retval = IDASetUserData(solver_, model_);
-            checkOutput(retval, "IDASetUserData");
+        // Set pointer to model data
+        retval = IDASetUserData(solver_, model_);
+        checkOutput(retval, "IDASetUserData");
 
-            // Set tolerances
-            realtype rtol;
-            realtype atol;
+        // Set tolerances
+        realtype rtol;
+        realtype atol;
 
-            model_->setTolerances(rtol, atol); ///< \todo Function name should be "getTolerances"!
-            retval = IDASStolerances(solver_, rtol, atol);
-            checkOutput(retval, "IDASStolerances");
+        model_->setTolerances(rtol, atol); ///< \todo Function name should be "getTolerances"!
+        retval = IDASStolerances(solver_, rtol, atol);
+        checkOutput(retval, "IDASStolerances");
+        
+        IdxT msa;
+        model_->setMaxSteps(msa);
 
-            IdxT msa;
-            model_->setMaxSteps(msa);
+        /// \todo Need to set max number of steps based on user input!
+        retval = IDASetMaxNumSteps(solver_, msa);
+        checkOutput(retval, "IDASetMaxNumSteps");
 
-            /// \todo Need to set max number of steps based on user input!
-            retval = IDASetMaxNumSteps(solver_, msa);
-            checkOutput(retval, "IDASetMaxNumSteps");
+        // Tag differential variables
+        std::vector<bool>& tag = model_->tag();
+        if (static_cast<IdxT>(tag.size()) == model_->size())
+        {
+            tag_ = N_VClone(yy_);
+            checkAllocation((void*) tag_, "N_VClone");
+            model_->tagDifferentiable();
+            copyVec(tag, tag_);
 
-            // Tag differential variables
-            std::vector<bool> &tag = model_->tag();
-            if (static_cast<IdxT>(tag.size()) == model_->size())
-            {
-                tag_ = N_VClone(yy_);
-                checkAllocation((void *)tag_, "N_VClone");
-                model_->tagDifferentiable();
-                copyVec(tag, tag_);
-
-                retval = IDASetId(solver_, tag_);
-                checkOutput(retval, "IDASetId");
-                retval = IDASetSuppressAlg(solver_, SUNTRUE);
-                checkOutput(retval, "IDASetSuppressAlg");
-            }
+            retval = IDASetId(solver_, tag_);
+            checkOutput(retval, "IDASetId");
+            retval = IDASetSuppressAlg(solver_, SUNTRUE);
+            checkOutput(retval, "IDASetSuppressAlg");
+        }
 
-            // Set up linear solver
-            this->configureLinearSolver();
+        // Set up linear solver
+        this->configureLinearSolver();
 
-            return retval;
-        }
+        return retval;
+    }
 
-        template <class ScalarT, typename IdxT>
-        int Ida<ScalarT, IdxT>::configureLinearSolver()
+    template <class ScalarT, typename IdxT>
+    int Ida<ScalarT, IdxT>::configureLinearSolver()
+    {
+        int retval = 0;
+        if (model_->hasJacobian())
         {
-            int retval = 0;
-            //Setup KLU for when Jacobian is needed
-            if (model_->hasJacobian())
-            {
-                JacobianMat_ = SUNSparseMatrix(model_->size(), model_->size(), model_->size() * model_->size(), CSR_MAT, context_);
-                checkAllocation((void *)JacobianMat_, "SUNSparseMatrix");
+            JacobianMat_ = SUNSparseMatrix(model_->size(), model_->size(), model_->size() * model_->size(), CSR_MAT, context_);
+            checkAllocation((void*) JacobianMat_, "SUNSparseMatrix");
 
-                linearSolver_ = SUNLinSol_KLU(yy_, JacobianMat_, context_);
-                checkAllocation((void *)linearSolver_, "SUNLinSol_KLU");
+            linearSolver_ = SUNLinSol_KLU(yy_, JacobianMat_, context_);
+            checkAllocation((void*) linearSolver_, "SUNLinSol_KLU");
 
-                retval = IDASetLinearSolver(solver_, linearSolver_, JacobianMat_);
-                checkOutput(retval, "IDASetLinearSolver");
+            retval = IDASetLinearSolver(solver_, linearSolver_, JacobianMat_);
+            checkOutput(retval, "IDASetLinearSolver");
 
-                retval = IDASetJacFn(solver_, this->Jac);
-                checkOutput(retval, "IDASetJacFn");
-            }
-            else
-            {
-                JacobianMat_ = SUNDenseMatrix(model_->size(), model_->size(), context_);
-                checkAllocation((void *)JacobianMat_, "SUNDenseMatrix");
+            retval = IDASetJacFn(solver_, this->Jac);
+            checkOutput(retval, "IDASetJacFn");
+        }
+        else
+        {
+            JacobianMat_ = SUNDenseMatrix(model_->size(), model_->size(), context_);
+            checkAllocation((void*) JacobianMat_, "SUNDenseMatrix");
 
-                linearSolver_ = SUNLinSol_Dense(yy_, JacobianMat_, context_);
-                checkAllocation((void *)linearSolver_, "SUNLinSol_Dense");
+            linearSolver_ = SUNLinSol_Dense(yy_, JacobianMat_, context_);
+            checkAllocation((void*) linearSolver_, "SUNLinSol_Dense");
 
-                retval = IDASetLinearSolver(solver_, linearSolver_, JacobianMat_);
-                checkOutput(retval, "IDASetLinearSolver");
-            }
+            retval = IDASetLinearSolver(solver_, linearSolver_, JacobianMat_);
+            checkOutput(retval, "IDASetLinearSolver");
 
-            return retval;
         }
 
-        template <class ScalarT, typename IdxT>
-        int Ida<ScalarT, IdxT>::getDefaultInitialCondition()
-        {
-            model_->initialize();
+        return retval;
+    }
 
-            copyVec(model_->y(), yy_);
-            copyVec(model_->yp(), yp_);
+    template <class ScalarT, typename IdxT>
+    int Ida<ScalarT, IdxT>::getDefaultInitialCondition()
+    {
+        model_->initialize();
 
-            return 0;
-        }
+        copyVec(model_->y(), yy_);
+        copyVec(model_->yp(), yp_);
 
-        template <class ScalarT, typename IdxT>
-        int Ida<ScalarT, IdxT>::setIntegrationTime(real_type t_init, real_type t_final, int nout)
-        {
-            t_init_ = t_init;
-            t_final_ = t_final;
-            nout_ = nout;
-            return 0;
-        }
+        return 0;
+    }
 
-        template <class ScalarT, typename IdxT>
-        int Ida<ScalarT, IdxT>::initializeSimulation(real_type t0, bool findConsistent)
-        {
-            int retval = 0;
+    template <class ScalarT, typename IdxT>
+    int Ida<ScalarT, IdxT>::setIntegrationTime(real_type t_init, real_type t_final, int nout)
+    {
+        t_init_  = t_init;
+        t_final_ = t_final;
+        nout_    = nout;
+        return 0;
+    }
+
+    template <class ScalarT, typename IdxT>
+    int Ida<ScalarT, IdxT>::initializeSimulation(real_type t0, bool findConsistent)
+    {
+        int retval = 0;
 
-            // Need to reinitialize IDA to set to get correct initial conditions
-            retval = IDAReInit(solver_, t0, yy_, yp_);
-            checkOutput(retval, "IDAReInit");
+        // Need to reinitialize IDA to set to get correct initial conditions
+        retval = IDAReInit(solver_, t0, yy_, yp_);
+        checkOutput(retval, "IDAReInit");
 
-            // Find a consistent set of initial conditions for DAE
-            if (findConsistent)
-            {
-                int initType = IDA_Y_INIT;
-
-                if (tag_)
-                    initType = IDA_YA_YDP_INIT;
+        // Find a consistent set of initial conditions for DAE
+        if (findConsistent)
+        {
+            int initType = IDA_Y_INIT;
 
-                retval = IDACalcIC(solver_, initType, 0.1);
-                checkOutput(retval, "IDACalcIC");
-            }
+            if (tag_)
+                initType = IDA_YA_YDP_INIT;
 
-            return retval;
+            retval = IDACalcIC(solver_, initType, 0.1);
+            checkOutput(retval, "IDACalcIC");
         }
 
-        template <class ScalarT, typename IdxT>
-        int Ida<ScalarT, IdxT>::runSimulation(real_type tf, int nout)
-        {
-            int retval = 0;
-            int iout = 0;
-            real_type tret;
-            real_type dt = tf / nout;
-            real_type tout = dt;
-
-            /* In loop, call IDASolve, print results, and test for error.
-             *     Break out of loop when NOUT preset output times have been reached. */
-            // printOutput(0.0);
-            while (nout > iout)
+        return retval;
+    }
+
+    template <class ScalarT, typename IdxT>
+    int Ida<ScalarT, IdxT>::runSimulation(real_type tf, int nout)
+    {
+        int retval = 0;
+        int iout = 0;
+        real_type tret;
+        real_type dt = tf/nout;
+        real_type tout = dt;
+
+        /* In loop, call IDASolve, print results, and test for error.
+         *     Break out of loop when NOUT preset output times have been reached. */
+        //printOutput(0.0);
+        while(nout > iout)
+        {
+            retval = IDASolve(solver_, tout, &tret, yy_, yp_, IDA_NORMAL);
+            checkOutput(retval, "IDASolve");
+            //printOutput(tout); 
+
+            if (retval == IDA_SUCCESS)
             {
-                retval = IDASolve(solver_, tout, &tret, yy_, yp_, IDA_NORMAL);
-                checkOutput(retval, "IDASolve");
-                // printOutput(tout);
-
-                if (retval == IDA_SUCCESS)
-                {
-                    ++iout;
-                    tout += dt;
-                }
+                ++iout;
+                tout += dt;
             }
-            // std::cout << "\n";
-            return retval;
         }
+        //std::cout << "\n";
+        return retval;
+    }
 
-        template <class ScalarT, typename IdxT>
-        int Ida<ScalarT, IdxT>::deleteSimulation()
-        {
-            IDAFree(&solver_);
-            SUNLinSolFree(linearSolver_);
-            N_VDestroy(yy_);
-            N_VDestroy(yp_);
-            return 0;
-        }
+    template <class ScalarT, typename IdxT>
+    int Ida<ScalarT, IdxT>::deleteSimulation()
+    {
+        IDAFree(&solver_);
+        SUNLinSolFree(linearSolver_);
+        N_VDestroy(yy_);
+        N_VDestroy(yp_);
+        return 0;
+    }
 
-        template <class ScalarT, typename IdxT>
-        int Ida<ScalarT, IdxT>::configureQuadrature()
-        {
-            int retval = 0;
 
-            // Create and initialize quadratures
-            q_ = N_VNew_Serial(model_->size_quad(), context_);
-            checkAllocation((void *)q_, "N_VNew_Serial");
+    template <class ScalarT, typename IdxT>
+    int Ida<ScalarT, IdxT>::configureQuadrature()
+    {
+        int retval = 0;
 
-            // Set integrand function and allocate quadrature workspace
-            retval = IDAQuadInit(solver_, this->Integrand, q_);
-            checkOutput(retval, "IDAQuadInit");
+        // Create and initialize quadratures
+        q_ = N_VNew_Serial(model_->size_quad(), context_);
+        checkAllocation((void*) q_, "N_VNew_Serial");
 
-            // Set tolerances and error control for quadratures
-            real_type rtol, atol;
-            model_->setTolerances(rtol, atol);
+        // Set integrand function and allocate quadrature workspace
+        retval = IDAQuadInit(solver_, this->Integrand, q_);
+        checkOutput(retval, "IDAQuadInit");
 
-            // Set tolerances for quadrature stricter than for integration
-            retval = IDAQuadSStolerances(solver_, rtol * 0.1, atol * 0.1);
-            checkOutput(retval, "IDAQuadSStolerances");
+        // Set tolerances and error control for quadratures
+        real_type rtol, atol;
+        model_->setTolerances(rtol, atol);
 
-            // Include quadrature in eror checking
-            retval = IDASetQuadErrCon(solver_, SUNTRUE);
-            checkOutput(retval, "IDASetQuadErrCon");
+        // Set tolerances for quadrature stricter than for integration
+        retval = IDAQuadSStolerances(solver_, rtol*0.1, atol*0.1);
+        checkOutput(retval, "IDAQuadSStolerances");
 
-            return retval;
-        }
+        // Include quadrature in eror checking
+        retval = IDASetQuadErrCon(solver_, SUNTRUE);
+        checkOutput(retval, "IDASetQuadErrCon");
 
-        template <class ScalarT, typename IdxT>
-        int Ida<ScalarT, IdxT>::initializeQuadrature()
-        {
-            int retval = 0;
+        return retval;
+    }
 
-            // Set all quadratures to zero
-            N_VConst(RCONST(0.0), q_);
 
-            // Initialize quadratures
-            retval = IDAQuadReInit(solver_, q_);
-            checkOutput(retval, "IDAQuadInit");
+    template <class ScalarT, typename IdxT>
+    int Ida<ScalarT, IdxT>::initializeQuadrature()
+    {
+        int retval = 0;
 
-            return retval;
-        }
+        // Set all quadratures to zero
+        N_VConst(RCONST(0.0), q_);
 
-        template <class ScalarT, typename IdxT>
-        int Ida<ScalarT, IdxT>::runSimulationQuadrature(real_type tf, int nout)
-        {
-            int retval = 0;
-            real_type tret;
+        // Initialize quadratures
+        retval = IDAQuadReInit(solver_, q_);
+        checkOutput(retval, "IDAQuadInit");
 
-            // std::cout << "Forward integration for initial value problem ... \n";
+        return retval;
+    }
 
-            real_type dt = tf / nout;
-            real_type tout = dt;
-            // printOutput(0.0);
-            // printSpecial(0.0, yy_);
-            for (int i = 0; i < nout; ++i)
+
+    template <class ScalarT, typename IdxT>
+    int Ida<ScalarT, IdxT>::runSimulationQuadrature(real_type tf, int nout)
+    {
+        int retval = 0;
+        real_type tret;
+
+        //std::cout << "Forward integration for initial value problem ... \n";
+
+        real_type dt = tf/nout;
+        real_type tout = dt;
+        //printOutput(0.0);
+        //printSpecial(0.0, yy_);
+        for(int i = 0; i < nout; ++i)
+        {
+            retval = IDASolve(solver_, tout, &tret, yy_, yp_, IDA_NORMAL);
+            checkOutput(retval, "IDASolve");
+            //printSpecial(tout, yy_);
+            //printOutput(tout); 
+
+            if (retval == IDA_SUCCESS)
             {
-                retval = IDASolve(solver_, tout, &tret, yy_, yp_, IDA_NORMAL);
-                checkOutput(retval, "IDASolve");
-                // printSpecial(tout, yy_);
-                // printOutput(tout);
-
-                if (retval == IDA_SUCCESS)
-                {
-                    tout += dt;
-                }
-
-                retval = IDAGetQuad(solver_, &tret, q_);
-                checkOutput(retval, "IDAGetQuad");
+                tout += dt;
             }
 
-            return retval;
+            retval = IDAGetQuad(solver_, &tret, q_);
+            checkOutput(retval, "IDAGetQuad");
         }
 
-        template <class ScalarT, typename IdxT>
-        int Ida<ScalarT, IdxT>::deleteQuadrature()
-        {
-            IDAQuadFree(solver_);
-            N_VDestroy(q_);
+        return retval;
+    }
 
-            return 0;
-        }
 
-        template <class ScalarT, typename IdxT>
-        int Ida<ScalarT, IdxT>::configureAdjoint()
-        {
-            // Allocate adjoint vector, derivatives and quadrature
-            yyB_ = N_VNew_Serial(model_->size(), context_);
-            checkAllocation((void *)yyB_, "N_VNew_Serial");
+    template <class ScalarT, typename IdxT>
+    int Ida<ScalarT, IdxT>::deleteQuadrature()
+    {
+        IDAQuadFree(solver_);
+        N_VDestroy(q_);
 
-            ypB_ = N_VClone(yyB_);
-            checkAllocation((void *)ypB_, "N_VClone");
+        return 0;
+    }
 
-            qB_ = N_VNew_Serial(model_->size_opt(), context_);
-            checkAllocation((void *)qB_, "N_VNew_Serial");
 
-            return 0;
-        }
+    template <class ScalarT, typename IdxT>
+    int Ida<ScalarT, IdxT>::configureAdjoint()
+    {
+        // Allocate adjoint vector, derivatives and quadrature
+      yyB_ = N_VNew_Serial(model_->size(), context_);
+        checkAllocation((void*) yyB_, "N_VNew_Serial");
 
-        template <class ScalarT, typename IdxT>
-        int Ida<ScalarT, IdxT>::initializeAdjoint(IdxT steps)
-        {
-            int retval = 0;
+        ypB_ = N_VClone(yyB_);
+        checkAllocation((void*) ypB_, "N_VClone");
 
-            // Create adjoint workspace
-            retval = IDAAdjInit(solver_, steps, IDA_HERMITE);
-            checkOutput(retval, "IDAAdjInit");
+        qB_ = N_VNew_Serial(model_->size_opt(), context_);
+        checkAllocation((void*) qB_, "N_VNew_Serial");
 
-            return retval;
-        }
+        return 0;
+    }
 
-        template <class ScalarT, typename IdxT>
-        int Ida<ScalarT, IdxT>::initializeBackwardSimulation(real_type tf)
-        {
-            int retval = 0;
-            realtype rtol;
-            realtype atol;
+    template <class ScalarT, typename IdxT>
+    int Ida<ScalarT, IdxT>::initializeAdjoint(IdxT steps)
+    {
+        int retval = 0;
+
+        // Create adjoint workspace
+        retval = IDAAdjInit(solver_, steps, IDA_HERMITE);
+        checkOutput(retval, "IDAAdjInit");
+
+        return retval;
+    }
+
+    template <class ScalarT, typename IdxT>
+    int Ida<ScalarT, IdxT>::initializeBackwardSimulation(real_type tf)
+    {
+        int retval = 0;
+        realtype rtol;
+        realtype atol;
 
-            model_->initializeAdjoint();
+        model_->initializeAdjoint();
 
-            copyVec(model_->yB(), yyB_);
-            copyVec(model_->ypB(), ypB_);
-            N_VConst(0.0, qB_);
+        copyVec(model_->yB(),  yyB_);
+        copyVec(model_->ypB(), ypB_);
+        N_VConst(0.0, qB_);
 
-            retval = IDACreateB(solver_, &backwardID_);
-            checkOutput(retval, "IDACreateB");
+        retval = IDACreateB(solver_, &backwardID_);
+        checkOutput(retval, "IDACreateB");
 
-            // IDAInitB must be called after forward simulation run.
-            retval = IDAInitB(solver_, backwardID_, this->adjointResidual, tf, yyB_, ypB_);
-            checkOutput(retval, "IDAInitB");
+        // IDAInitB must be called after forward simulation run.
+        retval = IDAInitB(solver_, backwardID_, this->adjointResidual, tf, yyB_, ypB_);
+        checkOutput(retval, "IDAInitB");
 
-            model_->setTolerances(rtol, atol);
-            retval = IDASStolerancesB(solver_, backwardID_, rtol, atol);
-            checkOutput(retval, "IDASStolerancesB");
+        model_->setTolerances(rtol, atol);
+        retval = IDASStolerancesB(solver_, backwardID_, rtol, atol);
+        checkOutput(retval, "IDASStolerancesB");
 
-            retval = IDASetUserDataB(solver_, backwardID_, model_);
-            checkOutput(retval, "IDASetUserDataB");
+        retval = IDASetUserDataB(solver_, backwardID_, model_);
+        checkOutput(retval, "IDASetUserDataB");
 
             /// \todo Need to set max number of steps based on user input!
             retval = IDASetMaxNumStepsB(solver_, backwardID_, 2000);
             checkOutput(retval, "IDASetMaxNumSteps");
 
-            // Set up linear solver
-            JacobianMatB_ = SUNDenseMatrix(model_->size(), model_->size(), context_);
-            checkAllocation((void *)JacobianMatB_, "SUNDenseMatrix");
+        // Set up linear solver
+        JacobianMatB_ = SUNDenseMatrix(model_->size(), model_->size(), context_);
+        checkAllocation((void*) JacobianMatB_, "SUNDenseMatrix");
 
-            linearSolverB_ = SUNLinSol_Dense(yyB_, JacobianMatB_, context_);
-            checkAllocation((void *)linearSolverB_, "SUNLinSol_Dense");
+        linearSolverB_ = SUNLinSol_Dense(yyB_, JacobianMatB_, context_);
+        checkAllocation((void*) linearSolverB_, "SUNLinSol_Dense");
 
-            retval = IDASetLinearSolverB(solver_, backwardID_, linearSolverB_, JacobianMatB_);
-            checkOutput(retval, "IDASetLinearSolverB");
+        retval = IDASetLinearSolverB(solver_, backwardID_, linearSolverB_, JacobianMatB_);
+        checkOutput(retval, "IDASetLinearSolverB");
 
-            // Also reinitialize quadratures.
-            retval = IDAQuadInitB(solver_, backwardID_, this->adjointIntegrand, qB_);
-            checkOutput(retval, "IDAQuadInitB");
 
-            // retval = IDAQuadSStolerancesB(solver_, backwardID_, rtol*1.1, atol*1.1);
-            retval = IDAQuadSStolerancesB(solver_, backwardID_, rtol * 0.1, atol * 0.1);
-            checkOutput(retval, "IDAQuadSStolerancesB");
+        // Also reinitialize quadratures.
+        retval = IDAQuadInitB(solver_, backwardID_, this->adjointIntegrand, qB_);
+        checkOutput(retval, "IDAQuadInitB");
 
-            // Include quadratures in error control
-            retval = IDASetQuadErrConB(solver_, backwardID_, SUNTRUE);
-            checkOutput(retval, "IDASetQuadErrConB");
+        //retval = IDAQuadSStolerancesB(solver_, backwardID_, rtol*1.1, atol*1.1);
+        retval = IDAQuadSStolerancesB(solver_, backwardID_, rtol*0.1, atol*0.1);
+        checkOutput(retval, "IDAQuadSStolerancesB");
 
-            return retval;
-        }
+        // Include quadratures in error control
+        retval = IDASetQuadErrConB(solver_, backwardID_, SUNTRUE);
+        checkOutput(retval, "IDASetQuadErrConB");
 
-        template <class ScalarT, typename IdxT>
-        int Ida<ScalarT, IdxT>::configureLinearSolverBackward()
-        {
-            int retval = 0;
 
-            // Create Jacobian matrix
-            JacobianMatB_ = SUNDenseMatrix(model_->size(), model_->size(), context_);
-            checkAllocation((void *)JacobianMatB_, "SUNDenseMatrix");
+        return retval;
+    }
 
-            // Create linear solver
-            linearSolverB_ = SUNLinSol_Dense(yyB_, JacobianMatB_, context_);
-            checkAllocation((void *)linearSolverB_, "SUNLinSol_Dense");
+    template <class ScalarT, typename IdxT>
+    int Ida<ScalarT, IdxT>::configureLinearSolverBackward()
+    {
+        int retval = 0;
 
-            // Attach linear solver to IDA
-            retval = IDASetLinearSolverB(solver_, backwardID_, linearSolverB_, JacobianMatB_);
-            checkOutput(retval, "IDASetLinearSolverB");
+        // Create Jacobian matrix
+        JacobianMatB_ = SUNDenseMatrix(model_->size(), model_->size(), context_);
+        checkAllocation((void*) JacobianMatB_, "SUNDenseMatrix");
 
-            return retval;
-        }
+        // Create linear solver
+        linearSolverB_ = SUNLinSol_Dense(yyB_, JacobianMatB_, context_);
+        checkAllocation((void*) linearSolverB_, "SUNLinSol_Dense");
 
-        template <class ScalarT, typename IdxT>
-        int Ida<ScalarT, IdxT>::runForwardSimulation(real_type tf, int nout)
-        {
-            int retval = 0;
-            int ncheck;
-            real_type time;
+        // Attach linear solver to IDA
+        retval = IDASetLinearSolverB(solver_, backwardID_, linearSolverB_, JacobianMatB_);
+        checkOutput(retval, "IDASetLinearSolverB");
 
-            // std::cout << "Forward integration for adjoint analysis ... \n";
+        return retval;
+    }
 
-            real_type dt = tf / nout;
-            real_type tout = dt;
-            for (int i = 0; i < nout; ++i)
-            {
-                retval = IDASolveF(solver_, tout, &time, yy_, yp_, IDA_NORMAL, &ncheck);
-                checkOutput(retval, "IDASolveF");
+    template <class ScalarT, typename IdxT>
+    int Ida<ScalarT, IdxT>::runForwardSimulation(real_type tf, int nout)
+    {
+        int retval = 0;
+        int ncheck;
+        real_type time;
+
+        //std::cout << "Forward integration for adjoint analysis ... \n";
 
-                if (retval == IDA_SUCCESS)
-                {
-                    tout += dt;
-                }
+        real_type dt = tf/nout;
+        real_type tout = dt;
+        for(int i = 0; i < nout; ++i)
+        {
+            retval = IDASolveF(solver_, tout, &time, yy_, yp_, IDA_NORMAL, &ncheck);
+            checkOutput(retval, "IDASolveF");
 
-                retval = IDAGetQuad(solver_, &time, q_);
-                checkOutput(retval, "IDASolve");
+            if (retval == IDA_SUCCESS)
+            {
+                tout += dt;
             }
 
-            return retval;
+            retval = IDAGetQuad(solver_, &time, q_);
+            checkOutput(retval, "IDASolve");
         }
 
-        template <class ScalarT, typename IdxT>
-        int Ida<ScalarT, IdxT>::runBackwardSimulation(real_type t_init)
-        {
-            int retval = 0;
-            long int nstB;
-            real_type time;
+        return retval;
+    }
 
-            // std::cout << "Backward integration for adjoint analysis ... ";
+    template <class ScalarT, typename IdxT>
+    int Ida<ScalarT, IdxT>::runBackwardSimulation(real_type t_init)
+    {
+        int retval = 0;
+        long int nstB;
+        real_type time;
 
-            retval = IDASolveB(solver_, t_init, IDA_NORMAL);
-            checkOutput(retval, "IDASolveB");
+        //std::cout << "Backward integration for adjoint analysis ... ";
 
-            IDAGetNumSteps(IDAGetAdjIDABmem(solver_, backwardID_), &nstB);
-            // std::cout << "done ( nst = " << nstB << " )\n";
+        retval = IDASolveB(solver_, t_init, IDA_NORMAL);
+        checkOutput(retval, "IDASolveB");
 
-            retval = IDAGetB(solver_, backwardID_, &time, yyB_, ypB_);
-            checkOutput(retval, "IDAGetB");
+        IDAGetNumSteps(IDAGetAdjIDABmem(solver_, backwardID_), &nstB);
+        //std::cout << "done ( nst = " << nstB << " )\n";
 
-            retval = IDAGetQuadB(solver_, backwardID_, &time, qB_);
-            checkOutput(retval, "IDAGetQuadB");
+        retval = IDAGetB(solver_, backwardID_, &time, yyB_, ypB_);
+        checkOutput(retval, "IDAGetB");
 
-            return retval;
-        }
+        retval = IDAGetQuadB(solver_, backwardID_, &time, qB_);
+        checkOutput(retval, "IDAGetQuadB");
 
-        template <class ScalarT, typename IdxT>
-        int Ida<ScalarT, IdxT>::deleteAdjoint()
-        {
-            IDAAdjFree(solver_);
-            return 0;
-        }
-
-        template <class ScalarT, typename IdxT>
-        int Ida<ScalarT, IdxT>::Residual(realtype tres, N_Vector yy, N_Vector yp, N_Vector rr, void *user_data)
-        {
-            ModelLib::ModelEvaluator<ScalarT, IdxT> *model = static_cast<ModelLib::ModelEvaluator<ScalarT, IdxT> *>(user_data);
+        return retval;
+    }
 
-            model->updateTime(tres, 0.0);
-            copyVec(yy, model->y());
-            copyVec(yp, model->yp());
+    template <class ScalarT, typename IdxT>
+    int Ida<ScalarT, IdxT>::deleteAdjoint()
+    {
+        IDAAdjFree(solver_);
+        return 0;
+    }
 
-            model->evaluateResidual();
-            const std::vector<ScalarT> &f = model->getResidual();
-            copyVec(f, rr);
+    template <class ScalarT, typename IdxT>
+    int Ida<ScalarT, IdxT>::Residual(realtype tres, N_Vector yy, N_Vector yp, N_Vector rr, void *user_data)
+    {
+        ModelLib::ModelEvaluator<ScalarT, IdxT>* model = static_cast<ModelLib::ModelEvaluator<ScalarT, IdxT>*>(user_data);
 
-            return 0;
-        }
+        model->updateTime(tres, 0.0);
+        copyVec(yy, model->y());
+        copyVec(yp, model->yp());
 
-        template <class ScalarT, typename IdxT>
-        int Ida<ScalarT, IdxT>::Jac(realtype t, realtype cj, N_Vector yy, N_Vector yp, N_Vector resvec, SUNMatrix J, void *user_data, N_Vector tmp1, N_Vector tmp2, N_Vector tmp3)
-        {
+        model->evaluateResidual();
+        const std::vector<ScalarT>& f = model->getResidual();
+        copyVec(f, rr);
 
-            ModelLib::ModelEvaluator<ScalarT, IdxT> *model = static_cast<ModelLib::ModelEvaluator<ScalarT, IdxT> *>(user_data);
+        return 0;
+    }
 
-            model->updateTime(t, cj);
-            copyVec(yy, model->y());
-            copyVec(yp, model->yp());
+    template <class ScalarT, typename IdxT>
+    int Ida<ScalarT, IdxT>::Jac(realtype t, realtype cj, N_Vector yy, N_Vector yp, N_Vector resvec, SUNMatrix J, void *user_data, N_Vector tmp1, N_Vector tmp2, N_Vector tmp3)
+    {
 
-            model->evaluateJacobian();
-            COO_Matrix<ScalarT, IdxT> Jac = model->getJacobian();
+        ModelLib::ModelEvaluator<ScalarT, IdxT>* model = static_cast<ModelLib::ModelEvaluator<ScalarT, IdxT>*>(user_data);
 
-            // Get reference to the jacobian entries
-            std::tuple<std::vector<IdxT> &, std::vector<IdxT> &, std::vector<ScalarT> &> tpm = Jac.getEntries();
-            const auto [r, c, val] = tpm;
 
-            // get the CSR row pointers from COO matrix
-            std::vector<IdxT> csrrowdata = Jac.getCSRRowData();
+        model->updateTime(t, cj);
+        copyVec(yy, model->y());
+        copyVec(yp, model->yp());
 
-            SUNMatZero(J);
+        model->evaluateJacobian();
+        COO_Matrix<ScalarT, IdxT> Jac = model->getJacobian();
+        
+        //Get reference to the jacobian entries
+        std::tuple<std::vector<IdxT>&, std::vector<IdxT>&, std::vector<ScalarT>&> tpm = Jac.getEntries();
+        const auto [r, c, val] = tpm;
 
-            // Set row pointers
-            sunindextype *rowptrs = SUNSparseMatrix_IndexPointers(J);
-            for (unsigned int i = 0; i < csrrowdata.size(); i++)
-            {
-                rowptrs[i] = csrrowdata[i];
-            }
+        //get the CSR row pointers from COO matrix
+        std::vector<IdxT> csrrowdata = Jac.getCSRRowData();
 
-            sunindextype *colvals = SUNSparseMatrix_IndexValues(J);
-            realtype *data = SUNSparseMatrix_Data(J);
-            // Copy data from model jac to sundials
-            for (unsigned int i = 0; i < c.size(); i++)
-            {
-                colvals[i] = c[i];
-                data[i] = val[i];
-            }
+        SUNMatZero(J);
 
-            return 0;
+        //Set row pointers
+        sunindextype *rowptrs = SUNSparseMatrix_IndexPointers(J);
+        for (unsigned int i = 0; i < csrrowdata.size() ; i++)
+        {
+            rowptrs[i] = csrrowdata[i];
         }
 
-        template <class ScalarT, typename IdxT>
-        int Ida<ScalarT, IdxT>::Integrand(realtype tt, N_Vector yy, N_Vector yp, N_Vector rhsQ, void *user_data)
+        sunindextype *colvals = SUNSparseMatrix_IndexValues(J);
+        realtype *data = SUNSparseMatrix_Data(J);
+        //Copy data from model jac to sundials
+        for (unsigned int i = 0; i < c.size(); i++ )
         {
-            ModelLib::ModelEvaluator<ScalarT, IdxT> *model = static_cast<ModelLib::ModelEvaluator<ScalarT, IdxT> *>(user_data);
+            colvals[i] = c[i];
+            data[i] = val[i];
+        }
 
-            model->updateTime(tt, 0.0);
-            copyVec(yy, model->y());
-            copyVec(yp, model->yp());
+        return 0;
+    }
 
-            model->evaluateIntegrand();
-            const std::vector<ScalarT> &g = model->getIntegrand();
-            copyVec(g, rhsQ);
+    template <class ScalarT, typename IdxT>
+    int Ida<ScalarT, IdxT>::Integrand(realtype tt, N_Vector yy, N_Vector yp, N_Vector rhsQ, void *user_data)
+    {
+        ModelLib::ModelEvaluator<ScalarT, IdxT>* model = static_cast<ModelLib::ModelEvaluator<ScalarT, IdxT>*>(user_data);
 
-            return 0;
-        }
+        model->updateTime(tt, 0.0);
+        copyVec(yy, model->y());
+        copyVec(yp, model->yp());
 
-        template <class ScalarT, typename IdxT>
-        int Ida<ScalarT, IdxT>::adjointResidual(realtype tt, N_Vector yy, N_Vector yp, N_Vector yyB, N_Vector ypB, N_Vector rrB, void *user_data)
-        {
-            ModelLib::ModelEvaluator<ScalarT, IdxT> *model = static_cast<ModelLib::ModelEvaluator<ScalarT, IdxT> *>(user_data);
+        model->evaluateIntegrand();
+        const std::vector<ScalarT>& g = model->getIntegrand();
+        copyVec(g, rhsQ);
 
-            model->updateTime(tt, 0.0);
-            copyVec(yy, model->y());
-            copyVec(yp, model->yp());
-            copyVec(yyB, model->yB());
-            copyVec(ypB, model->ypB());
+        return 0;
+    }
 
-            model->evaluateAdjointResidual();
-            const std::vector<ScalarT> &fB = model->getAdjointResidual();
-            copyVec(fB, rrB);
+    template <class ScalarT, typename IdxT>
+    int Ida<ScalarT, IdxT>::adjointResidual(realtype tt, N_Vector yy, N_Vector yp, N_Vector yyB, N_Vector ypB, N_Vector rrB, void *user_data)
+    {
+        ModelLib::ModelEvaluator<ScalarT, IdxT>* model = static_cast<ModelLib::ModelEvaluator<ScalarT, IdxT>*>(user_data);
 
-            return 0;
-        }
+        model->updateTime(tt, 0.0);
+        copyVec(yy, model->y());
+        copyVec(yp, model->yp());
+        copyVec(yyB, model->yB());
+        copyVec(ypB, model->ypB());
 
-        template <class ScalarT, typename IdxT>
-        int Ida<ScalarT, IdxT>::adjointIntegrand(realtype tt, N_Vector yy, N_Vector yp, N_Vector yyB, N_Vector ypB, N_Vector rhsQB, void *user_data)
-        {
-            ModelLib::ModelEvaluator<ScalarT, IdxT> *model = static_cast<ModelLib::ModelEvaluator<ScalarT, IdxT> *>(user_data);
+        model->evaluateAdjointResidual();
+        const std::vector<ScalarT>& fB = model->getAdjointResidual();
+        copyVec(fB, rrB);
 
-            model->updateTime(tt, 0.0);
-            copyVec(yy, model->y());
-            copyVec(yp, model->yp());
-            copyVec(yyB, model->yB());
-            copyVec(ypB, model->ypB());
+        return 0;
+    }
 
-            model->evaluateAdjointIntegrand();
-            const std::vector<ScalarT> &gB = model->getAdjointIntegrand();
-            copyVec(gB, rhsQB);
 
-            return 0;
-        }
+    template <class ScalarT, typename IdxT>
+    int Ida<ScalarT, IdxT>::adjointIntegrand(realtype tt, N_Vector yy, N_Vector yp, N_Vector yyB, N_Vector ypB, N_Vector rhsQB, void *user_data)
+    {
+        ModelLib::ModelEvaluator<ScalarT, IdxT>* model = static_cast<ModelLib::ModelEvaluator<ScalarT, IdxT>*>(user_data);
+
+        model->updateTime(tt, 0.0);
+        copyVec(yy, model->y());
+        copyVec(yp, model->yp());
+        copyVec(yyB, model->yB());
+        copyVec(ypB, model->ypB());
+
+        model->evaluateAdjointIntegrand();
+        const std::vector<ScalarT>& gB = model->getAdjointIntegrand();
+        copyVec(gB, rhsQB);
+
+        return 0;
+    }
+
 
-        template <class ScalarT, typename IdxT>
-        void Ida<ScalarT, IdxT>::copyVec(const N_Vector x, std::vector<ScalarT> &y)
+    template <class ScalarT, typename IdxT>
+    void Ida<ScalarT, IdxT>::copyVec(const N_Vector x, std::vector< ScalarT >& y)
+    {
+        const ScalarT* xdata = NV_DATA_S(x);
+        for(unsigned int i = 0; i < y.size(); ++i)
         {
-            const ScalarT *xdata = NV_DATA_S(x);
-            for (unsigned int i = 0; i < y.size(); ++i)
-            {
-                y[i] = xdata[i];
-            }
+            y[i] = xdata[i];
         }
+    }
+
 
-        template <class ScalarT, typename IdxT>
-        void Ida<ScalarT, IdxT>::copyVec(const std::vector<ScalarT> &x, N_Vector y)
+    template <class ScalarT, typename IdxT>
+    void Ida<ScalarT, IdxT>::copyVec(const std::vector< ScalarT >& x, N_Vector y)
+    {
+        ScalarT* ydata = NV_DATA_S(y);
+        for(unsigned int i = 0; i < x.size(); ++i)
         {
-            ScalarT *ydata = NV_DATA_S(y);
-            for (unsigned int i = 0; i < x.size(); ++i)
-            {
-                ydata[i] = x[i];
-            }
+            ydata[i] = x[i];
         }
+    }
 
-        template <class ScalarT, typename IdxT>
-        void Ida<ScalarT, IdxT>::copyVec(const std::vector<bool> &x, N_Vector y)
+    template <class ScalarT, typename IdxT>
+    void Ida<ScalarT, IdxT>::copyVec(const std::vector< bool >& x, N_Vector y)
+    {
+        ScalarT* ydata = NV_DATA_S(y);
+        for(unsigned int i = 0; i < x.size(); ++i)
         {
-            ScalarT *ydata = NV_DATA_S(y);
-            for (unsigned int i = 0; i < x.size(); ++i)
-            {
-                if (x[i])
-                    ydata[i] = 1.0;
-                else
-                    ydata[i] = 0.0;
-            }
+            if (x[i])
+                ydata[i] = 1.0;
+            else
+                ydata[i] = 0.0;
         }
+    }
 
-        template <class ScalarT, typename IdxT>
-        void Ida<ScalarT, IdxT>::printOutput(realtype t)
-        {
-            realtype *yval = N_VGetArrayPointer_Serial(yy_);
-            realtype *ypval = N_VGetArrayPointer_Serial(yp_);
 
-            std::cout << std::setprecision(5) << std::setw(7) << t << " ";
-            for (IdxT i = 0; i < model_->size(); ++i)
-            {
-                std::cout << yval[i] << " ";
-            }
-            for (IdxT i = 0; i < model_->size(); ++i)
-            {
-                std::cout << ypval[i] << " ";
-            }
-            std::cout << "\n";
-        }
+    template <class ScalarT, typename IdxT>
+    void Ida<ScalarT, IdxT>::printOutput(realtype t)
+    {
+        realtype *yval  = N_VGetArrayPointer_Serial(yy_);
+        realtype *ypval = N_VGetArrayPointer_Serial(yp_);
 
-        template <class ScalarT, typename IdxT>
-        void Ida<ScalarT, IdxT>::printSpecial(realtype t, N_Vector y)
+        std::cout << std::setprecision(5) << std::setw(7) << t << " ";
+        for (IdxT i = 0; i < model_->size(); ++i)
         {
-            realtype *yval = N_VGetArrayPointer_Serial(y);
-            IdxT N = static_cast<IdxT>(N_VGetLength_Serial(y));
-            std::cout << "{";
-            std::cout << std::setprecision(5) << std::setw(7) << t;
-            for (IdxT i = 0; i < N; ++i)
-            {
-                std::cout << ", " << yval[i];
-            }
-            std::cout << "},\n";
+            std::cout << yval[i] << " ";
+        }
+        for (IdxT i = 0; i < model_->size(); ++i)
+        {
+            std::cout << ypval[i] << " ";
         }
+        std::cout << "\n";
+    }
 
-        template <class ScalarT, typename IdxT>
-        void Ida<ScalarT, IdxT>::printFinalStats()
+    template <class ScalarT, typename IdxT>
+    void Ida<ScalarT, IdxT>::printSpecial(realtype t, N_Vector y)
+    {
+        realtype *yval = N_VGetArrayPointer_Serial(y);
+        IdxT N = static_cast<IdxT>(N_VGetLength_Serial(y));
+        std::cout << "{";
+        std::cout << std::setprecision(5) << std::setw(7) << t;
+        for (IdxT i = 0; i < N; ++i)
         {
-            int retval = 0;
-            void *mem = solver_;
-            long int nst;
-            long int nre;
-            long int nje;
-            long int nni;
-            long int netf;
-            long int ncfn;
-
-            retval = IDAGetNumSteps(mem, &nst);
-            checkOutput(retval, "IDAGetNumSteps");
-            retval = IDAGetNumResEvals(mem, &nre);
-            checkOutput(retval, "IDAGetNumResEvals");
-            retval = IDAGetNumJacEvals(mem, &nje);
-            checkOutput(retval, "IDAGetNumJacEvals");
-            retval = IDAGetNumNonlinSolvIters(mem, &nni);
-            checkOutput(retval, "IDAGetNumNonlinSolvIters");
-            retval = IDAGetNumErrTestFails(mem, &netf);
-            checkOutput(retval, "IDAGetNumErrTestFails");
-            retval = IDAGetNumNonlinSolvConvFails(mem, &ncfn);
-            checkOutput(retval, "IDAGetNumNonlinSolvConvFails");
-
-            // std::cout << "\nFinal Run Statistics: \n\n";
-            std::cout << "Number of steps                    = " << nst << "\n";
-            std::cout << "Number of residual evaluations     = " << nre << "\n";
-            // std::cout << "Number of Jacobian evaluations     = " <<  nje  << "\n";
-            std::cout << "Number of nonlinear iterations     = " << nni << "\n";
-            std::cout << "Number of error test failures      = " << netf << "\n";
-            std::cout << "Number of nonlinear conv. failures = " << ncfn << "\n";
+            std::cout << ", " << yval[i];
         }
+        std::cout << "},\n";
+    }
 
-        template <class ScalarT, typename IdxT>
-        void Ida<ScalarT, IdxT>::checkAllocation(void *v, const char *functionName)
+    template <class ScalarT, typename IdxT>
+    void Ida<ScalarT, IdxT>::printFinalStats()
+    {
+        int retval = 0;
+        void* mem = solver_;
+        long int nst;
+        long int nre;
+        long int nje;
+        long int nni;
+        long int netf;
+        long int ncfn;
+
+        retval = IDAGetNumSteps(mem, &nst);
+        checkOutput(retval, "IDAGetNumSteps");
+        retval = IDAGetNumResEvals(mem, &nre);
+        checkOutput(retval, "IDAGetNumResEvals");
+        retval = IDAGetNumJacEvals(mem, &nje);
+        checkOutput(retval, "IDAGetNumJacEvals");
+        retval = IDAGetNumNonlinSolvIters(mem, &nni);
+        checkOutput(retval, "IDAGetNumNonlinSolvIters");
+        retval = IDAGetNumErrTestFails(mem, &netf);
+        checkOutput(retval, "IDAGetNumErrTestFails");
+        retval = IDAGetNumNonlinSolvConvFails(mem, &ncfn);
+        checkOutput(retval, "IDAGetNumNonlinSolvConvFails");
+
+        // std::cout << "\nFinal Run Statistics: \n\n";
+        std::cout << "Number of steps                    = " <<  nst  << "\n";
+        std::cout << "Number of residual evaluations     = " <<  nre  << "\n";
+        //std::cout << "Number of Jacobian evaluations     = " <<  nje  << "\n";
+        std::cout << "Number of nonlinear iterations     = " <<  nni  << "\n";
+        std::cout << "Number of error test failures      = " <<  netf << "\n";
+        std::cout << "Number of nonlinear conv. failures = " <<  ncfn << "\n";
+    }
+
+
+    template <class ScalarT, typename IdxT>
+    void Ida<ScalarT, IdxT>::checkAllocation(void* v, const char* functionName)
+    {
+        if (v == NULL)
         {
-            if (v == NULL)
-            {
-                std::cerr << "\nERROR: Function " << functionName << " failed -- returned NULL pointer!\n\n";
-                throw SundialsException();
-            }
+            std::cerr << "\nERROR: Function " << functionName << " failed -- returned NULL pointer!\n\n";
+            throw SundialsException();
         }
+    }
 
-        template <class ScalarT, typename IdxT>
-        void Ida<ScalarT, IdxT>::checkOutput(int retval, const char *functionName)
+    template <class ScalarT, typename IdxT>
+    void Ida<ScalarT, IdxT>::checkOutput(int retval, const char* functionName)
+    {
+        if (retval < 0)
         {
-            if (retval < 0)
-            {
-                std::cerr << "\nERROR: Function " << functionName << " failed with flag " << retval << "!\n\n";
-                throw SundialsException();
-            }
+            std::cerr << "\nERROR: Function " << functionName << " failed with flag " << retval << "!\n\n";
+            throw SundialsException();
         }
+    }
 
-        // Compiler will prevent building modules with data type incompatible with realtype
-        template class Ida<realtype, long int>;
-        template class Ida<realtype, int>;
-        template class Ida<realtype, size_t>;
+    // Compiler will prevent building modules with data type incompatible with realtype
+    template class Ida<realtype, long int>;
+    template class Ida<realtype, int>;
+    template class Ida<realtype, size_t>;
 
-    } // namespace Sundials
+} // namespace Sundials
 } // namespace AnalysisManager