diff --git a/trick_source/trick_utils/SAIntegrator/examples/CannonBall/CannonBall.cpp b/trick_source/trick_utils/SAIntegrator/examples/CannonBall/CannonBall.cpp
index 0d00754e..05a0d80a 100644
--- a/trick_source/trick_utils/SAIntegrator/examples/CannonBall/CannonBall.cpp
+++ b/trick_source/trick_utils/SAIntegrator/examples/CannonBall/CannonBall.cpp
@@ -27,15 +27,9 @@ void calc_derivs( double t, double state[], double derivs[], void* udata) {
 double impact( double t, double state[], RootFinder* root_finder, void* udata) {
     double root_error = root_finder->find_roots(t, state[1]);
     if (root_error == 0.0) {
-        printf("---------------------------------------------------------------\n");
-        printf("Impact at t = %5.10f x = %5.10f y = %5.10f.\n", t, state[0], state[1]);
-        printf("---------------------------------------------------------------\n");
         root_finder->init();
-        state[1] = 0.0000000001; // pos_y
-        state[2] = 0.0;          // vel_x
-        state[3] = 0.0;          // vel_y
-        state[4] = 0.0;          // acc_x
-        state[5] = 0.0;          // acc_y
+        state[2] =  0.9 * state[2];
+        state[3] = -0.9 * state[3];
     }
     return (root_error);
 }
@@ -55,17 +49,14 @@ int main ( int argc, char* argv[]) {
     double t = 0.0;
     SA::RK2Integrator integ(dt, 6, state_var_p, state_var_p, calc_derivs, NULL);
 
-// Uncomment the following two lines to have the RootFinder detect impact.
-//  RootFinder root_finder(0.00000000001, Negative);
-//  integ.add_Rootfinder(&root_finder, &impact);
+    RootFinder root_finder(0.00000000001, Negative);
+    integ.add_Rootfinder(&root_finder, &impact);
 
     init_state(cannon);
     print_header();
     print_state( t, cannon);
-    while (t < 5.3) {
-        integ.load();
-        integ.step();
-        integ.unload();
+    while (t < 20.0) {
+        integ.integrate();
         t = integ.getIndyVar();
         print_state( t, cannon);
     }
diff --git a/trick_source/trick_utils/SAIntegrator/examples/DoubleIntegral/DoubleIntegral.cpp b/trick_source/trick_utils/SAIntegrator/examples/DoubleIntegral/DoubleIntegral.cpp
index 91822b07..c712177e 100644
--- a/trick_source/trick_utils/SAIntegrator/examples/DoubleIntegral/DoubleIntegral.cpp
+++ b/trick_source/trick_utils/SAIntegrator/examples/DoubleIntegral/DoubleIntegral.cpp
@@ -23,9 +23,7 @@ void deriv_x( double x, double state_x[], double derivs_x[], void* udata) {
     integ_y->set_out_vars(state_y);
     integ_y->setIndyVar(context->start);
     while (integ_y->getIndyVar() <= context->end) {
-        integ_y->load();
-        integ_y->step();
-        integ_y->unload();
+        integ_y->integrate();
     }
     derivs_x[0] = area; // volume = ∫ area dx
 }
@@ -40,9 +38,7 @@ double doubleIntegral( double x_start, double x_end, double y_start, double y_en
     SA::RK2Integrator integ_x (dx, 1, state_x, state_x, deriv_x, &y_integ_context);
     integ_x.setIndyVar(x_start);
     while (integ_x.getIndyVar()<=x_end) {
-        integ_x.load();
-        integ_x.step();
-        integ_x.unload();
+        integ_x.integrate();
     }
     return volume;
 }
diff --git a/trick_source/trick_utils/SAIntegrator/include/SAIntegrator.hh b/trick_source/trick_utils/SAIntegrator/include/SAIntegrator.hh
index 357dea95..30153ab1 100644
--- a/trick_source/trick_utils/SAIntegrator/include/SAIntegrator.hh
+++ b/trick_source/trick_utils/SAIntegrator/include/SAIntegrator.hh
@@ -7,12 +7,15 @@ namespace SA {
         double X_out;      // Independent Variable Out
         double default_h;  // Default step-size
         void*  user_data;  // User data
+        void advanceIndyVar();
     public:
         Integrator(double h, void* udata): X_in(0.0), X_out(0.0), default_h(h), user_data(udata) {};
         virtual ~Integrator() {}
-        virtual void step();
         virtual void load();
         virtual void unload();
+        virtual void step();
+        void integrate();
+        virtual double undo_integrate();
         double getIndyVar();
         void setIndyVar(double v);
     };
@@ -27,30 +30,47 @@ namespace SA {
         double** inVars;         // Array of pointers to the variables whose values comprise the input state vector.
         double** outVars;        // Array of pointers to the variables to which the output state vector values are to be disseminated.
         DerivsFunc derivs_func;  // Pointer to the function that calculates the derivatives to be integrated.
-        double last_h;           // What was the last step-size? For undo_step().
-
+     // void advanceIndyVar();   // Inherited from Integrator::advanceIndyVar()
     public:
         FirstOrderODEIntegrator( double h, int N, double* in_vars[], double* out_vars[], DerivsFunc dfunc, void* udata);
         virtual ~FirstOrderODEIntegrator();
-        virtual double undo_step();
-        void load();
-        void unload();
-        void load_from_outState();
-        double** set_in_vars( double* in_vars[]);
-        double** set_out_vars( double* out_vars[]);
+        void load();                               // Overrides Integrator::load(), Final
+        void unload();                             // Overrides Integrator::unload(), Final
+        virtual void step();                       // Overrides Integrator::step(). Virtual
+        // integrate();                               Inherited from Integrator::integrate(), Final
+        virtual double undo_integrate();           // Overrides Integrator::undo_integrate(), Virtual
+        void load_from_outState();                 // New, Final
+        double** set_in_vars( double* in_vars[]);  // New, Final
+        double** set_out_vars( double* out_vars[]);// New, Final
+        // double getIndyVar();
+        // void setIndyVar(double v);
     };
 
     typedef double (*RootErrorFunc)( double x, double state[], RootFinder* root_finder, void* udata);
 
     class FirstOrderODEVariableStepIntegrator : public FirstOrderODEIntegrator {
+    private:
         RootFinder* root_finder;
         RootErrorFunc root_error_func;
+    protected:
+        double last_h;
+        // void advanceIndyVar();                                           // Inherited from Integrator::advanceIndyVar(), Final
+        void advanceIndyVar( double h);
     public:
         FirstOrderODEVariableStepIntegrator( double h, int N, double* in_vars[], double* out_vars[], DerivsFunc dfunc, void* udata);
         ~FirstOrderODEVariableStepIntegrator();
-        virtual void variable_step( double h);
-        void add_Rootfinder( RootFinder* root_finder, RootErrorFunc rfunc);
-        void step();
+        // load();                                                             Inherited from FirstOrderODEIntegrator, Final
+        // unload();                                                           Inherited from FirstOrderODEIntegrator, Final
+        void step();                                                        // Overrides FirstOrderODEIntegrator::step(), Final
+        // integrate();                                                     // Inherited from Integrator
+        double undo_integrate();                                            // Overrides FirstOrderODEIntegrator::undo_integrate(), Final
+        // load_from_outState();                                            // Inherited from FirstOrderODEIntegrator, Final
+        // set_in_vars();                                                   // Inherited from FirstOrderODEIntegrator, Final
+        // set_out_vars();                                                  // Inherited from FirstOrderODEIntegrator, Final
+        virtual void variable_step( double h);                              // New, Should be overridden in derived classes, Virtual
+        void add_Rootfinder( RootFinder* root_finder, RootErrorFunc rfunc); // New, Final
+        // double getIndyVar();                                             // Inherited from Integrator
+        // void setIndyVar(double v);                                       // Inherited from Integrator
     private:
         void find_roots(double h, unsigned int depth);
     };
@@ -60,7 +80,19 @@ namespace SA {
         double *derivs;
         EulerIntegrator( double h, int N, double* in_vars[], double* out_vars[], DerivsFunc dfunc, void* udata);
         ~EulerIntegrator();
-        void variable_step( double h);
+        // load();                                                              Inherited from FirstOrderODEIntegrator
+        // unload();                                                            Inherited from FirstOrderODEIntegrator
+        // void step();                                                         Inherited from FirstOrderODEVariableStepIntegrator
+        // void integrate()                                                     Inherited from FirstOrderODEVariableStepIntegrator
+        // double undo_integrate();                                             Inherited from FirstOrderODEVariableStepIntegrator
+        // load_from_outState();                                                Inherited from FirstOrderODEIntegrator
+        // set_in_vars();                                                       Inherited from FirstOrderODEIntegrator
+        // set_out_vars();                                                      Inherited from FirstOrderODEIntegrator
+        void variable_step( double h);                                       // Overrides FirstOrderODEVariableStepIntegrator::variable_step()
+        // void integrate(double h)                                             Inherited from FirstOrderODEVariableStepIntegrator
+        // void add_Rootfinder( RootFinder* root_finder, RootErrorFunc rfunc);  Inherited from FirstOrderODEVariableStepIntegrator
+        // double getIndyVar();                                                 Inherited from Integrator
+        // void setIndyVar(double v);                                           Inherited from Integrator
     };
 
     class HeunsMethod : public FirstOrderODEVariableStepIntegrator {
@@ -69,7 +101,7 @@ namespace SA {
         double *derivs[2];
         HeunsMethod( double h, int N, double* in_vars[], double* out_vars[], DerivsFunc dfunc, void* udata);
         ~HeunsMethod();
-        void variable_step( double h);
+        void variable_step( double h); // Overrides FirstOrderODEVariableStepIntegrator::variable_step()
     };
 
     class RK2Integrator : public FirstOrderODEVariableStepIntegrator {
@@ -78,7 +110,7 @@ namespace SA {
         double *derivs[2];
         RK2Integrator( double h, int N, double* in_vars[], double* out_vars[], DerivsFunc derivs_func, void* udata);
         ~RK2Integrator();
-        void variable_step( double h);
+        void variable_step( double h); // Overrides FirstOrderODEVariableStepIntegrator::variable_step()
     };
 
     class RK4Integrator : public FirstOrderODEVariableStepIntegrator {
@@ -87,7 +119,7 @@ namespace SA {
         double *derivs[4];
         RK4Integrator( double h, int N, double* in_vars[], double* out_vars[], DerivsFunc dfunc, void* udata);
         ~RK4Integrator();
-        void variable_step( double h);
+        void variable_step( double h); // Overrides FirstOrderODEVariableStepIntegrator::variable_step()
     };
 
     class RK3_8Integrator : public FirstOrderODEVariableStepIntegrator {
@@ -96,18 +128,18 @@ namespace SA {
         double *derivs[4];
         RK3_8Integrator( double h, int N, double* in_vars[], double* out_vars[], DerivsFunc dfunc, void* udata);
         ~RK3_8Integrator();
-        void variable_step( double h);
+        void variable_step( double h); // Overrides FirstOrderODEVariableStepIntegrator::variable_step()
     };
 
     class EulerCromerIntegrator : public Integrator {
     protected:
-        unsigned int state_size;
-        double** x_p;
-        double** v_p;
-        double* x_in;
-        double* v_in;
-        double* x_out;
-        double* v_out;
+        unsigned int nDimensions;
+        double** pos_p;
+        double** vel_p;
+        double* pos_in;
+        double* vel_in;
+        double* pos_out;
+        double* vel_out;
         double* g_out;
         double* f_out;
         DerivsFunc gderivs;
@@ -120,25 +152,26 @@ namespace SA {
         void step();
         void load();
         void unload();
+        double undo_integrate();
     };
 
-    // class ABMIntegrator : public FirstOrderODEIntegrator {
-    // protected:
-    //     double **    deriv_history;
-    //     double *     composite_deriv;
-    //     unsigned int bufferSize;
-    //     unsigned int algorithmHistorySize;
-    //     unsigned int hix;
-    //     unsigned int currentHistorySize;
-    //
-    // public:
-    //     ABMIntegrator(int history_size, double h, int N, double* in_vars[], double* out_vars[], DerivsFunc dfunc, void* udata);
-    //     ~ABMIntegrator();
-    //     void step();
-    //     void undo_step();
-    // };
+    // AdamsBashforthMoulton Self Priming
+    class ABMIntegrator : public FirstOrderODEIntegrator {
+    protected:
+        double **    deriv_history;
+        double *     composite_deriv;
+        unsigned int bufferSize;
+        unsigned int algorithmHistorySize;
+        unsigned int hix;
+        unsigned int currentHistorySize;
+    public:
+        ABMIntegrator(int history_size, double h, int N, double* in_vars[], double* out_vars[], DerivsFunc dfunc, void* udata);
+        ~ABMIntegrator();
+        void step();
+        double undo_integrate();
+    };
 
-    // AdamsBashforthMoulton 2
+    // AdamsBashforthMoulton 2 with RK2 Priming
     class ABM2Integrator : public FirstOrderODEIntegrator {
     protected:
         double **    deriv_history;
@@ -147,15 +180,14 @@ namespace SA {
         unsigned int hix;
         unsigned int currentHistorySize;
         FirstOrderODEVariableStepIntegrator* priming_integrator;
-
     public:
         ABM2Integrator( double h, int N, double* in_vars[], double* out_vars[], DerivsFunc dfunc, void* udata);
         ~ABM2Integrator();
         void step();
-        double undo_step();
+        double undo_integrate();
     };
 
-    // AdamsBashforthMoulton 4
+    // AdamsBashforthMoulton 4 with RK4 Priming
     class ABM4Integrator : public FirstOrderODEIntegrator {
     protected:
         double **    deriv_history;
@@ -164,11 +196,10 @@ namespace SA {
         unsigned int hix;
         unsigned int currentHistorySize;
         FirstOrderODEVariableStepIntegrator* priming_integrator;
-
     public:
         ABM4Integrator( double h, int N, double* in_vars[], double* out_vars[], DerivsFunc dfunc, void* udata);
         ~ABM4Integrator();
         void step();
-        double undo_step();
+        double undo_integrate();
     };
 }
diff --git a/trick_source/trick_utils/SAIntegrator/src/SAIntegrator.cpp b/trick_source/trick_utils/SAIntegrator/src/SAIntegrator.cpp
index e52abc99..f93ecedb 100644
--- a/trick_source/trick_utils/SAIntegrator/src/SAIntegrator.cpp
+++ b/trick_source/trick_utils/SAIntegrator/src/SAIntegrator.cpp
@@ -3,11 +3,33 @@
 #include <iostream>
 #include "SAIntegrator.hh"
 
-void SA::Integrator::step() { X_out = X_in + default_h; }
-void SA::Integrator::load() { X_in = X_out; }
-void SA::Integrator::unload() {}
-double SA::Integrator::getIndyVar() { return X_out; }
-void SA::Integrator::setIndyVar(double v) { X_out = v; }
+void SA::Integrator::advanceIndyVar() {
+    X_out = X_in + default_h;
+}
+void SA::Integrator::step() {
+    advanceIndyVar();
+}
+void SA::Integrator::load() {
+    X_in = X_out;
+}
+void SA::Integrator::unload() {
+    std::cout << "SA::Integrator : Nothing to unload." << std::endl;
+}
+void SA::Integrator::integrate() {
+    load();
+    step();
+    unload();
+}
+double SA::Integrator::undo_integrate() {
+    X_out = X_in;
+    return (default_h);
+}
+double SA::Integrator::getIndyVar() {
+    return X_out;
+}
+void SA::Integrator::setIndyVar(double v) {
+    X_out = v;
+}
 
 // ------------------------------------------------------------
 SA::FirstOrderODEIntegrator::FirstOrderODEIntegrator(
@@ -59,13 +81,24 @@ double** SA::FirstOrderODEIntegrator::set_out_vars( double* out_vars[]){
     outVars = out_vars;
     return (ret);
 }
-double SA::FirstOrderODEIntegrator::undo_step() {
-        for (unsigned int i=0 ; i<state_size; i++ ) {
-            *(inVars[i]) = inState[i];
-        }
-        X_out = X_in;
-    return (last_h);
+double SA::FirstOrderODEIntegrator::undo_integrate() {
+    for (unsigned int i=0 ; i<state_size; i++ ) {
+        *(inVars[i]) = inState[i];
+    }
+    for (unsigned int i=0 ; i<state_size; i++ ) {
+        outState[i] = inState[i];
+    }
+    return (SA::Integrator::undo_integrate());
 }
+void SA::FirstOrderODEIntegrator::step() {
+    double derivs[state_size];
+    (*derivs_func)( X_in, inState, derivs, user_data);
+    for (unsigned int i=0; i<state_size; i++) {
+        outState[i] = inState[i] + derivs[i] * default_h;
+    }
+    advanceIndyVar();
+}
+
 // ------------------------------------------------------------
 SA::FirstOrderODEVariableStepIntegrator::FirstOrderODEVariableStepIntegrator(
     double h, int N, double* in_vars[], double* out_vars[], DerivsFunc dfunc, void* udata)
@@ -73,6 +106,10 @@ SA::FirstOrderODEVariableStepIntegrator::FirstOrderODEVariableStepIntegrator(
     root_finder = (RootFinder*) NULL;
     root_error_func = (RootErrorFunc)NULL;
 }
+double SA::FirstOrderODEVariableStepIntegrator::undo_integrate() {
+    SA::FirstOrderODEIntegrator::undo_integrate();
+    return (last_h);
+}
 SA::FirstOrderODEVariableStepIntegrator::~FirstOrderODEVariableStepIntegrator() {
     // User is responsible for deleting the root_finder.
 }
@@ -81,15 +118,23 @@ void SA::FirstOrderODEVariableStepIntegrator::add_Rootfinder(
     root_finder = rootFinder;
     root_error_func = rfunc;
 }
-void SA::FirstOrderODEVariableStepIntegrator::variable_step( double h) {
+void SA::FirstOrderODEVariableStepIntegrator::advanceIndyVar(double h) {
     last_h = h; X_out = X_in + h;
 }
+void SA::FirstOrderODEVariableStepIntegrator::variable_step( double h) {
+    double derivs[state_size];
+    (*derivs_func)( X_in, inState, derivs, user_data);
+    for (unsigned int i=0; i<state_size; i++) {
+        outState[i] = inState[i] + derivs[i] * h;
+    }
+    advanceIndyVar(h);
+}
 void SA::FirstOrderODEVariableStepIntegrator::find_roots(double h, unsigned int depth) {
     double new_h;
     double h_correction = (*root_error_func)( getIndyVar(), outState, root_finder, user_data );
     if (h_correction < 0.0) {
         while (h_correction != 0.0) {
-            new_h = undo_step() + h_correction;
+            new_h = undo_integrate() + h_correction;
             variable_step(new_h);
             h_correction = (*root_error_func)( getIndyVar(), outState, root_finder, user_data );
         }
@@ -106,6 +151,7 @@ void SA::FirstOrderODEVariableStepIntegrator::step() {
         find_roots( default_h, 5 );
     }
 }
+
 // ------------------------------------------------------------
 SA::EulerIntegrator::EulerIntegrator(
     double h, int N, double* in_vars[], double* out_vars[], DerivsFunc dfunc, void* udata)
@@ -120,7 +166,7 @@ void SA::EulerIntegrator::variable_step( double h) {
     for (unsigned int i=0; i<state_size; i++) {
         outState[i] = inState[i] + derivs[i] * h;
     }
-    SA::FirstOrderODEVariableStepIntegrator::variable_step(h);
+    advanceIndyVar(h);
 }
 // ------------------------------------------------------------
 SA::HeunsMethod::HeunsMethod(
@@ -146,7 +192,7 @@ void SA::HeunsMethod::variable_step( double h) {
     for (unsigned int i=0; i<state_size; i++) {
         outState[i] = inState[i] + (h/2) * ( derivs[0][i] +  derivs[1][i] );
     }
-    SA::FirstOrderODEVariableStepIntegrator::variable_step(h);
+    advanceIndyVar(h);
 }
 // ------------------------------------------------------------
 SA::RK2Integrator::RK2Integrator(
@@ -172,7 +218,7 @@ void SA::RK2Integrator::variable_step( double h) {
     for (unsigned int i=0; i<state_size; i++) {
         outState[i] = inState[i] + h * derivs[1][i];
     }
-    SA::FirstOrderODEVariableStepIntegrator::variable_step(h);
+    advanceIndyVar(h);
 }
 // ------------------------------------------------------------
 SA::RK4Integrator::RK4Integrator(
@@ -213,7 +259,7 @@ void SA::RK4Integrator::variable_step( double h) {
                                  (1/3.0)* derivs[2][i] +
                                  (1/6.0)* derivs[3][i]) * h;
     }
-    SA::FirstOrderODEVariableStepIntegrator::variable_step(h);
+    advanceIndyVar(h);
 }
 // ------------------------------------------------------------
 SA::RK3_8Integrator::RK3_8Integrator(
@@ -254,138 +300,146 @@ void SA::RK3_8Integrator::variable_step( double h) {
                                  (3/8.0)* derivs[2][i] +
                                  (1/8.0)* derivs[3][i]) * h;
     }
-    SA::FirstOrderODEVariableStepIntegrator::variable_step(h);
+    advanceIndyVar(h);
 }
 // ------------------------------------------------------------
 SA::EulerCromerIntegrator::EulerCromerIntegrator(
     double dt, int N, double* xp[], double* vp[], DerivsFunc gfunc, DerivsFunc ffunc, void* udata)
     :Integrator(dt, udata) {
-    state_size  = N;
-    x_p   = new double*[N];
-    v_p   = new double*[N];
-    x_in  = new double[N];
-    v_in  = new double[N];
-    g_out = new double[N];
-    f_out = new double[N];
-    x_out = new double[N];
-    v_out = new double[N];
+    nDimensions  = N;
+    pos_p   = new double*[N];
+    pos_in  = new double[N];
+    pos_out = new double[N];
+    vel_p   = new double*[N];
+    vel_in  = new double[N];
+    vel_out = new double[N];
+    g_out   = new double[N];
+    f_out   = new double[N];
+
     for (int i=0 ; i<N; i++ ) {
-        x_p[i] = xp[i];
-        v_p[i] = vp[i];
+        pos_p[i] = xp[i];
+        vel_p[i] = vp[i];
     }
     gderivs = gfunc;
     fderivs = ffunc;
 }
 SA::EulerCromerIntegrator::~EulerCromerIntegrator() {
-    delete x_p;
-    delete v_p;
-    delete x_in;
-    delete v_in;
+
+    delete pos_in;
+    delete pos_out;
+    delete vel_in;
+    delete vel_out;
     delete g_out;
     delete f_out;
-    delete x_out;
-    delete v_out;
 }
 void SA::EulerCromerIntegrator::step( double dt) {
-    (*gderivs)( X_in, x_in, g_out, user_data);
-    for (int i=0 ; i<state_size; i++ ) {
-        v_out[i] = v_in[i] + g_out[i] * dt;
+    (*gderivs)( X_in, pos_in, g_out, user_data);
+    for (int i=0 ; i<nDimensions; i++ ) {
+        vel_out[i] = vel_in[i] + g_out[i] * dt;
     }
-    (*fderivs)( X_in, v_in, f_out, user_data);
-    for (int i=0 ; i<state_size; i++ ) {
-        x_out[i] = x_in[i] + f_out[i] * dt;
+    (*fderivs)( X_in, vel_in, f_out, user_data);
+    for (int i=0 ; i<nDimensions; i++ ) {
+        pos_out[i] = pos_in[i] + f_out[i] * dt;
     }
     X_out = X_in + dt;
 }
 void SA::EulerCromerIntegrator::step() {
     step(default_h);
+    // if ((root_finder != NULL) && (root_error_func != NULL)) {
+    //     find_roots( default_h, 5 );
+    // }
 }
 void SA::EulerCromerIntegrator::load() {
-    for (int i=0 ; i<state_size; i++ ) {
-        x_in[i] = *(x_p[i]);
-        v_in[i] = *(v_p[i]);
+    for (int i=0 ; i<nDimensions; i++ ) {
+        pos_in[i] = *(pos_p[i]);
+        vel_in[i] = *(vel_p[i]);
     }
 }
 void SA::EulerCromerIntegrator::unload(){
-    for (int i=0 ; i<state_size; i++ ) {
-        *(x_p[i]) = x_out[i];
-        *(v_p[i]) = v_out[i];
+    for (int i=0 ; i<nDimensions; i++ ) {
+        *(pos_p[i]) = pos_out[i];
+        *(vel_p[i]) = vel_out[i];
     }
 }
+double SA::EulerCromerIntegrator::undo_integrate() {
+    for (int i=0 ; i<nDimensions; i++ ) {
+        *(pos_p[i]) = pos_in[i];
+        *(vel_p[i]) = vel_in[i];
+    }
+    return (SA::Integrator::undo_integrate());
+}
 
 // // Adams-Bashforth Coefficients
-// static const double ABCoeffs[5][5] = {
-//   {1.0,           0.0,           0.0,          0.0,           0.0},
-//   {(3/2.0),      (-1/2.0),       0.0,          0.0,           0.0},
-//   {(23/12.0),    (-16/12.0),    (5/12.0),      0.0,           0.0},
-//   {(55/24.0),    (-59/24.0),    (37/24.0),    (-9/24.0),      0.0},
-//   {(1901/720.0), (-2774/720.0), (2616/720.0), (-1274/720.0), (251/720.0)}
-// };
+static const double ABCoeffs[5][5] = {
+  {1.0,           0.0,           0.0,          0.0,           0.0},
+  {(3/2.0),      (-1/2.0),       0.0,          0.0,           0.0},
+  {(23/12.0),    (-16/12.0),    (5/12.0),      0.0,           0.0},
+  {(55/24.0),    (-59/24.0),    (37/24.0),    (-9/24.0),      0.0},
+  {(1901/720.0), (-2774/720.0), (2616/720.0), (-1274/720.0), (251/720.0)}
+};
 //
 // // Adams-Moulton Coefficients
-// static const double AMCoeffs[5][5] = {
-//   {1.0,          0.0,         0.0,          0.0,         0.0},
-//   {(1/2.0),     (1/2.0),      0.0,          0.0,         0.0},
-//   {(5/12.0),    (8/12.0),    (-1/12.0),     0.0,         0.0},
-//   {(9/24.0),    (19/24.0),   (-5/24.0),    (1/24.0),     0.0},
-//   {(251/720.0), (646/720.0), (-264/720.0), (106/720.0), (-19/720.0)}
-// };
-// SA::ABMIntegrator::ABMIntegrator ( int history_size, double h, int N, double* in_vars[], double* out_vars[], DerivsFunc func, void* udata)
-//     : FirstOrderODEIntegrator(h, N, in_vars, out_vars ,func, udata) {
-//
-//     algorithmHistorySize=history_size; // The amount of history that we intend to use in this ABM integrator.
-//     bufferSize=algorithmHistorySize+1; // The size of the buffer needs to be one more than the history so that an integration step can be reset.
-//     hix=0;
-//     currentHistorySize=0;            // How much derivative history is currently stored int the history buffer. Initially there will be none until we've primed the integrator.
-//     deriv_history = new double*[bufferSize];
-//     for (unsigned int i=0; i<bufferSize ; i++) {
-//         deriv_history[i] = new double[state_size];
-//     }
-//     composite_deriv = new double[state_size];
-// }
-// SA::ABMIntegrator::~ABMIntegrator() {
-//     for (int i=0; i<bufferSize ; i++) {
-//         delete (deriv_history[i]);
-//     }
-//     delete(deriv_history);
-//     delete(composite_deriv);
-// }
-// void SA::ABMIntegrator::step() {
-//
-//     hix = (hix+1)%bufferSize;                          // Move history index forward
-//     // (*derivs_func)( indyVar, inState, deriv_history[hix], user_data);
-//     (*derivs_func)( X_in, inState, deriv_history[hix], user_data); // Calculated and store the deriv for current, corrected state.
-//     // Increase the size of the stored history, up to the limit specified by the user.
-//     currentHistorySize = (currentHistorySize < algorithmHistorySize) ? currentHistorySize+1 : algorithmHistorySize;
-//     // (Predictor) Predict the next state using the Adams-Bashforth explicit method.
-//     for (int i=0; i<state_size; i++) {
-//         composite_deriv[i] = 0.0;
-//         for (int n=0,j=hix; n<currentHistorySize ; n++,j=(j+bufferSize-1)%bufferSize) {
-//             composite_deriv[i] += ABCoeffs[currentHistorySize-1][n] * deriv_history[j][i];
-//         }
-//         outState[i] = inState[i] + default_h * composite_deriv[i];
-//     }
-//     // Move history index forward, so we can temporarily store the derivative of the predicted next state.
-//     // We do not increase the currentHistorySize.
-//     hix = (hix+1)%bufferSize;
-//     (*derivs_func)( X_out, outState, deriv_history[hix], user_data); // Calc deriv for predicted next state.
-//     // (Corrector) Refine the next state using the Adams-Moulton implicit method. This is the corrected next state.
-//     for (int i=0; i<state_size; i++) {
-//         composite_deriv[i] = 0.0;
-//         for (int n=0,j=hix; n<currentHistorySize ; n++,j=(j+bufferSize-1)%bufferSize) {
-//             composite_deriv[i] += AMCoeffs[currentHistorySize-1][n] * deriv_history[j][i];
-//         }
-//         outState[i] = inState[i] + default_h * composite_deriv[i];
-//     }
-//     // Move history index backward, so we over-write the predicted state with the corrected state on our next step().
-//     hix = (hix+bufferSize-1)%bufferSize;
-//     SA::Integrator::step();
-// }
-//
-// double SA::ABMIntegrator::undo_step() {
-//     hix = (hix+bufferSize-1)%bufferSize;
-//     return (FirstOrderODEIntegrator::undo_step());
-// }
+static const double AMCoeffs[5][5] = {
+  {1.0,          0.0,         0.0,          0.0,         0.0},
+  {(1/2.0),     (1/2.0),      0.0,          0.0,         0.0},
+  {(5/12.0),    (8/12.0),    (-1/12.0),     0.0,         0.0},
+  {(9/24.0),    (19/24.0),   (-5/24.0),    (1/24.0),     0.0},
+  {(251/720.0), (646/720.0), (-264/720.0), (106/720.0), (-19/720.0)}
+};
+SA::ABMIntegrator::ABMIntegrator ( int history_size, double h, int N, double* in_vars[], double* out_vars[], DerivsFunc func, void* udata)
+    : FirstOrderODEIntegrator(h, N, in_vars, out_vars ,func, udata) {
+
+    algorithmHistorySize=history_size; // The amount of history that we intend to use in this ABM integrator.
+    bufferSize=algorithmHistorySize+1; // The size of the buffer needs to be one more than the history so that an integration step can be reset.
+    hix=0;
+    currentHistorySize=0;              // How much derivative history is currently stored int the history buffer. Initially there will be none until we've primed the integrator.
+    deriv_history = new double*[bufferSize];
+    for (unsigned int i=0; i<bufferSize ; i++) {
+        deriv_history[i] = new double[state_size];
+    }
+    composite_deriv = new double[state_size];
+}
+SA::ABMIntegrator::~ABMIntegrator() {
+    for (int i=0; i<bufferSize ; i++) {
+        delete (deriv_history[i]);
+    }
+    delete(deriv_history);
+    delete(composite_deriv);
+}
+void SA::ABMIntegrator::step() {
+    hix = (hix+1)%bufferSize;                                      // Move history index forward
+    (*derivs_func)( X_in, inState, deriv_history[hix], user_data); // Calculated and store the deriv for current, corrected state.
+    // Increase the size of the stored history, up to the limit specified by the user.
+    currentHistorySize = (currentHistorySize < algorithmHistorySize) ? currentHistorySize+1 : algorithmHistorySize;
+    // (Predictor) Predict the next state using the Adams-Bashforth explicit method.
+    for (int i=0; i<state_size; i++) {
+        composite_deriv[i] = 0.0;
+        for (int n=0,j=hix; n<currentHistorySize ; n++,j=(j+bufferSize-1)%bufferSize) {
+            composite_deriv[i] += ABCoeffs[currentHistorySize-1][n] * deriv_history[j][i];
+        }
+        outState[i] = inState[i] + default_h * composite_deriv[i];
+    }
+    // Move history index forward, so we can temporarily store the derivative of the predicted next state.
+    // We do not increase the currentHistorySize.
+    hix = (hix+1)%bufferSize;
+    (*derivs_func)( X_out, outState, deriv_history[hix], user_data); // Calc deriv for predicted next state.
+    // (Corrector) Refine the next state using the Adams-Moulton implicit method. This is the corrected next state.
+    for (int i=0; i<state_size; i++) {
+        composite_deriv[i] = 0.0;
+        for (int n=0,j=hix; n<currentHistorySize ; n++,j=(j+bufferSize-1)%bufferSize) {
+            composite_deriv[i] += AMCoeffs[currentHistorySize-1][n] * deriv_history[j][i];
+        }
+        outState[i] = inState[i] + default_h * composite_deriv[i];
+    }
+    // Move history index backward, so we over-write the predicted state with the corrected state on our next step().
+    hix = (hix+bufferSize-1)%bufferSize;
+    SA::Integrator::step();
+}
+
+double SA::ABMIntegrator::undo_integrate() {
+    hix = (hix+bufferSize-1)%bufferSize;
+    return (FirstOrderODEIntegrator::undo_integrate());
+}
 
 // =======================================================================
 static const double AB2Coeffs[2] = {(3/2.0), (-1/2.0)};
@@ -404,7 +458,7 @@ SA::ABM2Integrator::ABM2Integrator (
     }
     composite_deriv = new double[state_size];
 
-    // NEW: Create a priming integrator.
+    // Create a priming integrator.
     double** priming_integrator_in_vars = new double*[state_size];
     for (unsigned int i=0; i<N ; i++) {
         priming_integrator_in_vars[i] = &(inState[i]);
@@ -458,11 +512,11 @@ void SA::ABM2Integrator::step() {
         // Move history index backward, so we over-write the predicted state with the corrected state on our next step().
         hix = (hix+bufferSize-1)%bufferSize;
     }
-    SA::Integrator::step();
+    advanceIndyVar();
 }
-double SA::ABM2Integrator::undo_step() {
+double SA::ABM2Integrator::undo_integrate() {
     hix = (hix+bufferSize-1)%bufferSize;
-    return ( FirstOrderODEIntegrator::undo_step());
+    return ( FirstOrderODEIntegrator::undo_integrate());
 }
 // =======================================================================
 static const double AB4Coeffs[4] = {(55/24.0), (-59/24.0), (37/24.0), (-9/24.0)};
@@ -534,9 +588,9 @@ void SA::ABM4Integrator::step() {
         // Move history index backward, so we over-write the predicted state with the corrected state on our next step().
         hix = (hix+bufferSize-1)%bufferSize;
     }
-    SA::Integrator::step();
+    advanceIndyVar();
 }
-double SA::ABM4Integrator::undo_step() {
+double SA::ABM4Integrator::undo_integrate() {
     hix = (hix+bufferSize-1)%bufferSize;
-    return (FirstOrderODEIntegrator::undo_step());
+    return (FirstOrderODEIntegrator::undo_integrate());
 }
diff --git a/trick_source/trick_utils/SAIntegrator/unittest/FirstOrderODEIntegrator_unittest.cc b/trick_source/trick_utils/SAIntegrator/unittest/FirstOrderODEIntegrator_unittest.cc
new file mode 100644
index 00000000..022358f8
--- /dev/null
+++ b/trick_source/trick_utils/SAIntegrator/unittest/FirstOrderODEIntegrator_unittest.cc
@@ -0,0 +1,159 @@
+#include <gtest/gtest.h>
+#include <iostream>
+#include "SAIntegrator.hh"
+#include <math.h>
+
+#define EXCEPTABLE_ERROR 0.00000000001
+
+void calc_derivs( double x       __attribute__((unused)),
+                  double state[] __attribute__((unused)),
+                  double derivs[],
+                  void*  udata   __attribute__((unused))) {
+    derivs[0] = 1.0;
+    derivs[1] = 1.0;
+    derivs[2] = 1.0;
+    derivs[3] = 1.0;
+}
+
+/*
+The default integration algorithm for FirstOrderODEIntegrator is Euler integration.
+Tests:
+    * getIndyVar
+    * integrate
+        * load
+        * step
+        * unload
+*/
+TEST(FirstOrderODEIntegrator_unittest, test_1) {
+    double h = 0.1;
+    unsigned int N = 4;
+    double vars[4] = {1.0, 2.0, 3.0, 4.0} ;
+    double* varptrs[4] = {&vars[0], &vars[1], &vars[2], &vars[3]};
+    double x;
+
+    SA::FirstOrderODEIntegrator integ( h, N, varptrs, varptrs, calc_derivs, NULL);
+
+    x = integ.getIndyVar();
+    EXPECT_NEAR(x, 0.0, EXCEPTABLE_ERROR);
+
+    integ.integrate();
+    EXPECT_NEAR(vars[0], 1.1, EXCEPTABLE_ERROR);
+    EXPECT_NEAR(vars[1], 2.1, EXCEPTABLE_ERROR);
+    EXPECT_NEAR(vars[2], 3.1, EXCEPTABLE_ERROR);
+    EXPECT_NEAR(vars[3], 4.1, EXCEPTABLE_ERROR);
+
+    x = integ.getIndyVar();
+    EXPECT_NEAR(x, 0.1, EXCEPTABLE_ERROR);
+}
+/*
+Tests:
+    * getIndyVar
+    * integrate
+        * load
+        * step
+        * unload
+    * undo_integrate
+*/
+TEST(FirstOrderODEIntegrator_unittest, test_2) {
+    double h = 0.1;
+    unsigned int N = 4;
+    double vars[4] = {1.0, 2.0, 3.0, 4.0} ;
+    double* varptrs[4] = {&vars[0], &vars[1], &vars[2], &vars[3]};
+    double x;
+
+    SA::FirstOrderODEIntegrator integ( h, N, varptrs, varptrs, calc_derivs, NULL);
+
+    integ.integrate();
+    EXPECT_NEAR(vars[0], 1.1, EXCEPTABLE_ERROR);
+    EXPECT_NEAR(vars[1], 2.1, EXCEPTABLE_ERROR);
+    EXPECT_NEAR(vars[2], 3.1, EXCEPTABLE_ERROR);
+    EXPECT_NEAR(vars[3], 4.1, EXCEPTABLE_ERROR);
+
+    x = integ.getIndyVar();
+    EXPECT_NEAR(x, 0.1, EXCEPTABLE_ERROR);
+
+    integ.undo_integrate();
+    EXPECT_NEAR(vars[0], 1.0, EXCEPTABLE_ERROR);
+    EXPECT_NEAR(vars[1], 2.0, EXCEPTABLE_ERROR);
+    EXPECT_NEAR(vars[2], 3.0, EXCEPTABLE_ERROR);
+    EXPECT_NEAR(vars[3], 4.0, EXCEPTABLE_ERROR);
+
+    x = integ.getIndyVar();
+    EXPECT_NEAR(x, 0.0, EXCEPTABLE_ERROR);
+}
+
+/*
+Tests:
+    * getIndyVar
+    * integrate
+        * load
+        * step
+        * unload
+    * load_from_outState
+*/
+TEST(FirstOrderODEIntegrator_unittest, test_3) {
+    double h = 0.1;
+    unsigned int N = 4;
+    double vars[4] = {1.0, 2.0, 3.0, 4.0};
+    double* varptrs[4] = {&vars[0], &vars[1], &vars[2], &vars[3]};
+
+    SA::FirstOrderODEIntegrator integ( h, N, varptrs, varptrs, calc_derivs, NULL);
+
+    integ.integrate();
+    for (int i=0; i<5 ; i++) {
+        integ.load_from_outState();
+        integ.step();
+    }
+    integ.unload();
+    EXPECT_NEAR(vars[0], 1.6, EXCEPTABLE_ERROR);
+    EXPECT_NEAR(vars[1], 2.6, EXCEPTABLE_ERROR);
+    EXPECT_NEAR(vars[2], 3.6, EXCEPTABLE_ERROR);
+    EXPECT_NEAR(vars[3], 4.6, EXCEPTABLE_ERROR);
+
+    double x = integ.getIndyVar();
+    EXPECT_NEAR(x, 0.6, EXCEPTABLE_ERROR);
+}
+
+/*
+Tests:
+    * getIndyVar
+    * integrate
+        * load
+        * step
+        * unload
+    * load_from_outState
+    * set_out_vars
+*/
+TEST(FirstOrderODEIntegrator_unittest, test_4) {
+    double h = 0.1;
+    unsigned int N = 4;
+
+    double vars[4] = {1.0, 2.0, 3.0, 4.0};
+    double* varptrs[4] = {&vars[0], &vars[1], &vars[2], &vars[3]};
+    double altvars[4];
+    double* altvarptrs[4] = {&altvars[0], &altvars[1], &altvars[2], &altvars[3]};
+
+    SA::FirstOrderODEIntegrator integ( h, N, varptrs, varptrs, calc_derivs, NULL);
+
+    integ.integrate();
+    for (int i=0; i<5 ; i++) {
+        integ.load_from_outState();
+        integ.step();
+    }
+    integ.set_out_vars( altvarptrs);
+    integ.unload();
+    EXPECT_NEAR(altvars[0], 1.6, EXCEPTABLE_ERROR);
+    EXPECT_NEAR(altvars[1], 2.6, EXCEPTABLE_ERROR);
+    EXPECT_NEAR(altvars[2], 3.6, EXCEPTABLE_ERROR);
+    EXPECT_NEAR(altvars[3], 4.6, EXCEPTABLE_ERROR);
+
+    integ.set_out_vars( varptrs);
+    integ.unload();
+    EXPECT_NEAR(vars[0], 1.6, EXCEPTABLE_ERROR);
+    EXPECT_NEAR(vars[1], 2.6, EXCEPTABLE_ERROR);
+    EXPECT_NEAR(vars[2], 3.6, EXCEPTABLE_ERROR);
+    EXPECT_NEAR(vars[3], 4.6, EXCEPTABLE_ERROR);
+
+    double x = integ.getIndyVar();
+    EXPECT_NEAR(x, 0.6, EXCEPTABLE_ERROR);
+}
diff --git a/trick_source/trick_utils/SAIntegrator/unittest/FirstOrderODEVariableStepIntegrator_unittest.cc b/trick_source/trick_utils/SAIntegrator/unittest/FirstOrderODEVariableStepIntegrator_unittest.cc
new file mode 100644
index 00000000..14675eac
--- /dev/null
+++ b/trick_source/trick_utils/SAIntegrator/unittest/FirstOrderODEVariableStepIntegrator_unittest.cc
@@ -0,0 +1,172 @@
+#include <gtest/gtest.h>
+#include <iostream>
+#include "SAIntegrator.hh"
+#include <math.h>
+
+void deriv4( double t,
+             double state[] __attribute__((unused)),
+             double derivs[],
+             void* udata __attribute__((unused))) {
+    double t2 = t*t;
+    double t3 = t2*t;
+    derivs[0] = 3.0*t3 - 3.0*t2 - 3.0*t +  4.0;
+    derivs[1] = 2.0*t3 + 1.0*t2 - 5.0*t +  7.0;
+    derivs[2] =     t3 +     t2 + 5.0*t +  4.0;
+    derivs[3] =     t3 - 6.0*t2 + 4.0*t + 12.0;
+}
+
+void deriv2( double t,
+             double state[] __attribute__((unused)),
+             double derivs[],
+             void* udata __attribute__((unused))) {
+
+    derivs[0] = -2.0*t + 3.0;
+    derivs[1] =  5.0*t + 4.0;
+    derivs[2] =  3.0*t - 3.0;
+    derivs[3] =  5.0*t + 4.0;
+}
+
+void deriv1( double t __attribute__((unused)),
+             double state[] __attribute__((unused)),
+             double derivs[],
+             void* udata __attribute__((unused))) {
+    derivs[0] = 4.0;
+    derivs[1] = -4.0;
+    derivs[2] = M_PI;
+    derivs[3] = -M_PI;
+}
+
+#define EXCEPTABLE_ERROR 0.00000000001
+
+TEST(FirstOrderODEVariableStepIntegrator_unittest, EulerIntegrator_1) {
+
+   double state[4] = {0.0, 0.0, 0.0, 0.0};
+   double* state_var_p[4] = { &(state[0]), &(state[1]), &(state[2]), &(state[3])};
+   double dt = 0.01;
+   unsigned int count = 0;
+   double t = count * dt;
+   SA::EulerIntegrator integ(dt, 4, state_var_p, state_var_p, deriv1, NULL);
+   while (t < 2.0) {
+        integ.load();
+        integ.step();
+        integ.unload();
+        count ++;
+        t = count * dt;
+    }
+    EXPECT_NEAR(state[0], 8.0, EXCEPTABLE_ERROR);
+    EXPECT_NEAR(state[1], -8.0, EXCEPTABLE_ERROR);
+    EXPECT_NEAR(state[2], (2*M_PI), EXCEPTABLE_ERROR);
+    EXPECT_NEAR(state[3], (-2*M_PI), EXCEPTABLE_ERROR);
+}
+
+TEST(FirstOrderODEVariableStepIntegrator_unittest, RungeKutta2_1) {
+
+   double state[4] = {0.0, 0.0, 0.0, 0.0};
+   double* state_var_p[4] = { &(state[0]), &(state[1]), &(state[2]), &(state[3]) };
+   double dt = 0.01;
+   unsigned int count = 0;
+   double t = count * dt;
+   SA::RK2Integrator integ(dt, 4, state_var_p, state_var_p, deriv2, NULL);
+   while (t < 2.0) {
+        integ.load();
+        integ.step();
+        integ.unload();
+        count ++;
+        t = count * dt;
+    }
+    EXPECT_NEAR(state[0], 2.0, EXCEPTABLE_ERROR);
+    EXPECT_NEAR(state[1], 18.0, EXCEPTABLE_ERROR);
+    EXPECT_NEAR(state[2], 0.0, EXCEPTABLE_ERROR);
+    EXPECT_NEAR(state[3], 18.0, EXCEPTABLE_ERROR);
+}
+
+TEST(FirstOrderODEVariableStepIntegrator_unittest, RungeKutta4_1) {
+
+   double state[4] = {0.0, 0.0, 0.0, 0.0};
+   double* state_var_p[4] = { &(state[0]), &(state[1]), &(state[2]), &(state[3]) };
+   double dt = 0.01;
+   unsigned int count = 0;
+   double t = count * dt;
+   SA::RK4Integrator integ(dt, 4, state_var_p, state_var_p, deriv4, NULL);
+   while (t < 2.0) {
+        integ.load();
+        integ.step();
+        integ.unload();
+        count ++;
+        t = count * dt;
+    }
+    EXPECT_NEAR(state[0], 6.0, EXCEPTABLE_ERROR);
+    EXPECT_NEAR(state[1], (14.0+2/3.0), EXCEPTABLE_ERROR);
+    EXPECT_NEAR(state[2], (24+2/3.0), EXCEPTABLE_ERROR);
+    EXPECT_NEAR(state[3], 20.0, EXCEPTABLE_ERROR);
+}
+
+TEST(FirstOrderODEVariableStepIntegrator_unittest, RungeKutta38_1) {
+
+   double state[4] = {0.0, 0.0, 0.0, 0.0};
+   double* state_var_p[4] = { &(state[0]), &(state[1]), &(state[2]), &(state[3]) };
+   double dt = 0.01;
+   unsigned int count = 0;
+   double t = count * dt;
+   SA::RK3_8Integrator integ(dt, 4, state_var_p, state_var_p, deriv4, NULL);
+   while (t < 2.0) {
+        integ.load();
+        integ.step();
+        integ.unload();
+        count ++;
+        t = count * dt;
+    }
+    EXPECT_NEAR(state[0], 6.0, EXCEPTABLE_ERROR);
+    EXPECT_NEAR(state[1], (14.0+2/3.0), EXCEPTABLE_ERROR);
+    EXPECT_NEAR(state[2], (24+2/3.0), EXCEPTABLE_ERROR);
+    EXPECT_NEAR(state[3], 20.0, EXCEPTABLE_ERROR);
+}
+
+TEST(FirstOrderODEVariableStepIntegrator_unittest, EulerIntegrator_undo_integrate) {
+
+   double state[4] = {0.0, 0.0, 0.0, 0.0};
+   double* state_var_p[4] = { &(state[0]), &(state[1]), &(state[2]), &(state[3])};
+   double dx = 0.1;
+   double x;
+
+   SA::EulerIntegrator integ(dx, 4, state_var_p, state_var_p, deriv1, NULL);
+   integ.load();
+   integ.step();
+   integ.unload();
+
+   x = integ.getIndyVar();
+   EXPECT_NEAR(state[0],   (4.0/10.0), EXCEPTABLE_ERROR);
+   EXPECT_NEAR(state[1],  (-4.0/10.0), EXCEPTABLE_ERROR);
+   EXPECT_NEAR(state[2],  (M_PI/10.0), EXCEPTABLE_ERROR);
+   EXPECT_NEAR(state[3], (-M_PI/10.0), EXCEPTABLE_ERROR);
+   EXPECT_NEAR(x, 0.1, EXCEPTABLE_ERROR);
+
+   integ.undo_integrate();
+
+   x = integ.getIndyVar();
+   EXPECT_NEAR(state[0],  0.0, EXCEPTABLE_ERROR);
+   EXPECT_NEAR(state[1],  0.0, EXCEPTABLE_ERROR);
+   EXPECT_NEAR(state[2],  0.0, EXCEPTABLE_ERROR);
+   EXPECT_NEAR(state[3],  0.0, EXCEPTABLE_ERROR);
+   EXPECT_NEAR(x, 0.0, EXCEPTABLE_ERROR);
+
+   integ.undo_integrate();
+
+   x = integ.getIndyVar();
+   EXPECT_NEAR(state[0],  0.0, EXCEPTABLE_ERROR);
+   EXPECT_NEAR(state[1],  0.0, EXCEPTABLE_ERROR);
+   EXPECT_NEAR(state[2],  0.0, EXCEPTABLE_ERROR);
+   EXPECT_NEAR(state[3],  0.0, EXCEPTABLE_ERROR);
+   EXPECT_NEAR(x, 0.0, EXCEPTABLE_ERROR);
+
+   integ.load();
+   integ.step();
+   integ.unload();
+
+   x = integ.getIndyVar();
+   EXPECT_NEAR(state[0],   (4.0/10.0), EXCEPTABLE_ERROR);
+   EXPECT_NEAR(state[1],  (-4.0/10.0), EXCEPTABLE_ERROR);
+   EXPECT_NEAR(state[2],  (M_PI/10.0), EXCEPTABLE_ERROR);
+   EXPECT_NEAR(state[3], (-M_PI/10.0), EXCEPTABLE_ERROR);
+   EXPECT_NEAR(x, 0.1, EXCEPTABLE_ERROR);
+}
diff --git a/trick_source/trick_utils/SAIntegrator/unittest/Makefile b/trick_source/trick_utils/SAIntegrator/unittest/Makefile
index 5bf083dd..1a078d94 100644
--- a/trick_source/trick_utils/SAIntegrator/unittest/Makefile
+++ b/trick_source/trick_utils/SAIntegrator/unittest/Makefile
@@ -15,8 +15,10 @@ LIBDIRS = -L${SAI_LIBDIR} -L${GTEST_HOME}/lib64 -L${GTEST_HOME}/lib
 
 all: test
 
-test: SAIntegrator_unittest RootFinder_unittest
+test: SAIntegrator_unittest FirstOrderODEIntegrator_unittest FirstOrderODEVariableStepIntegrator_unittest RootFinder_unittest
 	./SAIntegrator_unittest --gtest_output=xml:${TRICK_HOME}/trick_test/SAIntegrator_unittest.xml
+	./FirstOrderODEIntegrator_unittest --gtest_output=xml:${TRICK_HOME}/trick_test/FirstOrderODEIntegrator_unittest.xml
+	./FirstOrderODEVariableStepIntegrator_unittest --gtest_output=xml:${TRICK_HOME}/trick_test/FirstOrderODEVariableStepIntegrator_unittest.xml
 	./RootFinder_unittest --gtest_output=xml:${TRICK_HOME}/trick_test/RootFinder_unittest.xml
 
 SAIntegrator_unittest.o : SAIntegrator_unittest.cc
@@ -25,6 +27,18 @@ SAIntegrator_unittest.o : SAIntegrator_unittest.cc
 SAIntegrator_unittest : ${SAI_LIBDIR}/${SAI_LIBNAME} SAIntegrator_unittest.o
 	$(TRICK_CXX) $(TRICK_CPPFLAGS) -o $@ $^ ${LIBDIRS} -lSAInteg -lgtest -lgtest_main -lpthread
 
+FirstOrderODEIntegrator_unittest.o : FirstOrderODEIntegrator_unittest.cc
+	$(TRICK_CXX) $(TRICK_CPPFLAGS) $(INCLUDE_DIRS) -c $<
+
+FirstOrderODEIntegrator_unittest : ${SAI_LIBDIR}/${SAI_LIBNAME} FirstOrderODEIntegrator_unittest.o
+	$(TRICK_CXX) $(TRICK_CPPFLAGS) -o $@ $^ ${LIBDIRS} -lSAInteg -lgtest -lgtest_main -lpthread
+
+FirstOrderODEVariableStepIntegrator_unittest.o : FirstOrderODEVariableStepIntegrator_unittest.cc
+	$(TRICK_CXX) $(TRICK_CPPFLAGS) $(INCLUDE_DIRS) -c $<
+
+FirstOrderODEVariableStepIntegrator_unittest : ${SAI_LIBDIR}/${SAI_LIBNAME} FirstOrderODEVariableStepIntegrator_unittest.o
+	$(TRICK_CXX) $(TRICK_CPPFLAGS) -o $@ $^ ${LIBDIRS} -lSAInteg -lgtest -lgtest_main -lpthread
+
 RootFinder_unittest.o : RootFinder_unittest.cc
 	$(TRICK_CXX) $(TRICK_CPPFLAGS) $(INCLUDE_DIRS) -c $<
 
@@ -38,4 +52,4 @@ clean:
 	${RM} *.o
 
 spotless: clean
-	${RM} SAIntegrator_unittest
+	${RM} FirstOrderODEVariableStepIntegrator_unittest
diff --git a/trick_source/trick_utils/SAIntegrator/unittest/SAIntegrator_unittest.cc b/trick_source/trick_utils/SAIntegrator/unittest/SAIntegrator_unittest.cc
index 7764ae2b..96c2a0d9 100644
--- a/trick_source/trick_utils/SAIntegrator/unittest/SAIntegrator_unittest.cc
+++ b/trick_source/trick_utils/SAIntegrator/unittest/SAIntegrator_unittest.cc
@@ -3,170 +3,72 @@
 #include "SAIntegrator.hh"
 #include <math.h>
 
-void deriv4( double t,
-             double state[] __attribute__((unused)),
-             double derivs[],
-             void* udata __attribute__((unused))) {
-    double t2 = t*t;
-    double t3 = t2*t;
-    derivs[0] = 3.0*t3 - 3.0*t2 - 3.0*t +  4.0;
-    derivs[1] = 2.0*t3 + 1.0*t2 - 5.0*t +  7.0;
-    derivs[2] =     t3 +     t2 + 5.0*t +  4.0;
-    derivs[3] =     t3 - 6.0*t2 + 4.0*t + 12.0;
-}
-
-void deriv2( double t,
-             double state[] __attribute__((unused)),
-             double derivs[],
-             void* udata __attribute__((unused))) {
-
-    derivs[0] = -2.0*t + 3.0;
-    derivs[1] =  5.0*t + 4.0;
-    derivs[2] =  3.0*t - 3.0;
-    derivs[3] =  5.0*t + 4.0;
-}
-
-void deriv1( double t __attribute__((unused)),
-             double state[] __attribute__((unused)),
-             double derivs[],
-             void* udata __attribute__((unused))) {
-    derivs[0] = 4.0;
-    derivs[1] = -4.0;
-    derivs[2] = M_PI;
-    derivs[3] = -M_PI;
-}
-
 #define EXCEPTABLE_ERROR 0.00000000001
 
-TEST(SAIntegrator_unittest, EulerIntegrator_1) {
+TEST(Integrator, test_getIndyVar_and_setIndyVar) {
+    double h = 0.1;
+    void* udata = NULL;
+    double x;
 
-   double state[4] = {0.0, 0.0, 0.0, 0.0};
-   double* state_var_p[4] = { &(state[0]), &(state[1]), &(state[2]), &(state[3])};
-   double dt = 0.01;
-   unsigned int count = 0;
-   double t = count * dt;
-   SA::EulerIntegrator integ(dt, 4, state_var_p, state_var_p, deriv1, NULL);
-   while (t < 2.0) {
-        integ.load();
-        integ.step();
-        integ.unload();
-        count ++;
-        t = count * dt;
-    }
-    EXPECT_NEAR(state[0], 8.0, EXCEPTABLE_ERROR);
-    EXPECT_NEAR(state[1], -8.0, EXCEPTABLE_ERROR);
-    EXPECT_NEAR(state[2], (2*M_PI), EXCEPTABLE_ERROR);
-    EXPECT_NEAR(state[3], (-2*M_PI), EXCEPTABLE_ERROR);
+    SA::Integrator  integ( h, udata);
+    //Confirm that initial value of the independent variable is 0.0;
+    x = integ.getIndyVar();
+    EXPECT_NEAR(x, 0.0, EXCEPTABLE_ERROR);
+
+    //Set the independent variable, and then confirm that it's the value we set.;
+    integ.setIndyVar(3.0);
+    x = integ.getIndyVar();
+    EXPECT_NEAR(x, 3.0, EXCEPTABLE_ERROR);
+}
+TEST(Integrator, test_load_and_step) {
+    double h = 0.1;
+    void* udata = NULL;
+    double x;
+
+    SA::Integrator  integ( h, udata);
+    integ.setIndyVar(4.0);
+
+    integ.load();
+    integ.step();
+    x = integ.getIndyVar();
+    EXPECT_NEAR(x, 4.1, EXCEPTABLE_ERROR);
+
+    integ.load();
+    integ.step();
+    x = integ.getIndyVar();
+    EXPECT_NEAR(x, 4.2, EXCEPTABLE_ERROR);
 }
 
-TEST(SAIntegrator_unittest, RungeKutta2_1) {
+TEST(Integrator, test_integrate) {
+    double h = 0.1;
+    void* udata = NULL;
+    double x;
 
-   double state[4] = {0.0, 0.0, 0.0, 0.0};
-   double* state_var_p[4] = { &(state[0]), &(state[1]), &(state[2]), &(state[3]) };
-   double dt = 0.01;
-   unsigned int count = 0;
-   double t = count * dt;
-   SA::RK2Integrator integ(dt, 4, state_var_p, state_var_p, deriv2, NULL);
-   while (t < 2.0) {
-        integ.load();
-        integ.step();
-        integ.unload();
-        count ++;
-        t = count * dt;
-    }
-    EXPECT_NEAR(state[0], 2.0, EXCEPTABLE_ERROR);
-    EXPECT_NEAR(state[1], 18.0, EXCEPTABLE_ERROR);
-    EXPECT_NEAR(state[2], 0.0, EXCEPTABLE_ERROR);
-    EXPECT_NEAR(state[3], 18.0, EXCEPTABLE_ERROR);
+    SA::Integrator  integ( h, udata);
+    integ.setIndyVar(4.0);
+
+    integ.integrate();
+    x = integ.getIndyVar();
+    EXPECT_NEAR(x, 4.1, EXCEPTABLE_ERROR);
+
+    integ.integrate();
+    x = integ.getIndyVar();
+    EXPECT_NEAR(x, 4.2, EXCEPTABLE_ERROR);
 }
 
-TEST(SAIntegrator_unittest, RungeKutta4_1) {
+TEST(Integrator, test_undo_integrate) {
+    double h = 0.1;
+    void* udata = NULL;
+    double x;
 
-   double state[4] = {0.0, 0.0, 0.0, 0.0};
-   double* state_var_p[4] = { &(state[0]), &(state[1]), &(state[2]), &(state[3]) };
-   double dt = 0.01;
-   unsigned int count = 0;
-   double t = count * dt;
-   SA::RK4Integrator integ(dt, 4, state_var_p, state_var_p, deriv4, NULL);
-   while (t < 2.0) {
-        integ.load();
-        integ.step();
-        integ.unload();
-        count ++;
-        t = count * dt;
-    }
-    EXPECT_NEAR(state[0], 6.0, EXCEPTABLE_ERROR);
-    EXPECT_NEAR(state[1], (14.0+2/3.0), EXCEPTABLE_ERROR);
-    EXPECT_NEAR(state[2], (24+2/3.0), EXCEPTABLE_ERROR);
-    EXPECT_NEAR(state[3], 20.0, EXCEPTABLE_ERROR);
-}
-
-TEST(SAIntegrator_unittest, RungeKutta38_1) {
-
-   double state[4] = {0.0, 0.0, 0.0, 0.0};
-   double* state_var_p[4] = { &(state[0]), &(state[1]), &(state[2]), &(state[3]) };
-   double dt = 0.01;
-   unsigned int count = 0;
-   double t = count * dt;
-   SA::RK3_8Integrator integ(dt, 4, state_var_p, state_var_p, deriv4, NULL);
-   while (t < 2.0) {
-        integ.load();
-        integ.step();
-        integ.unload();
-        count ++;
-        t = count * dt;
-    }
-    EXPECT_NEAR(state[0], 6.0, EXCEPTABLE_ERROR);
-    EXPECT_NEAR(state[1], (14.0+2/3.0), EXCEPTABLE_ERROR);
-    EXPECT_NEAR(state[2], (24+2/3.0), EXCEPTABLE_ERROR);
-    EXPECT_NEAR(state[3], 20.0, EXCEPTABLE_ERROR);
-}
-
-TEST(SAIntegrator_unittest, EulerIntegrator_undo_step) {
-
-   double state[4] = {0.0, 0.0, 0.0, 0.0};
-   double* state_var_p[4] = { &(state[0]), &(state[1]), &(state[2]), &(state[3])};
-   double dx = 0.1;
-   double x;
-
-   SA::EulerIntegrator integ(dx, 4, state_var_p, state_var_p, deriv1, NULL);
-   integ.load();
-   integ.step();
-   integ.unload();
-
-   x = integ.getIndyVar();
-   EXPECT_NEAR(state[0],   (4.0/10.0), EXCEPTABLE_ERROR);
-   EXPECT_NEAR(state[1],  (-4.0/10.0), EXCEPTABLE_ERROR);
-   EXPECT_NEAR(state[2],  (M_PI/10.0), EXCEPTABLE_ERROR);
-   EXPECT_NEAR(state[3], (-M_PI/10.0), EXCEPTABLE_ERROR);
-   EXPECT_NEAR(x, 0.1, EXCEPTABLE_ERROR);
-
-   integ.undo_step();
-
-   x = integ.getIndyVar();
-   EXPECT_NEAR(state[0],  0.0, EXCEPTABLE_ERROR);
-   EXPECT_NEAR(state[1],  0.0, EXCEPTABLE_ERROR);
-   EXPECT_NEAR(state[2],  0.0, EXCEPTABLE_ERROR);
-   EXPECT_NEAR(state[3],  0.0, EXCEPTABLE_ERROR);
-   EXPECT_NEAR(x, 0.0, EXCEPTABLE_ERROR);
-
-   integ.undo_step();
-
-   x = integ.getIndyVar();
-   EXPECT_NEAR(state[0],  0.0, EXCEPTABLE_ERROR);
-   EXPECT_NEAR(state[1],  0.0, EXCEPTABLE_ERROR);
-   EXPECT_NEAR(state[2],  0.0, EXCEPTABLE_ERROR);
-   EXPECT_NEAR(state[3],  0.0, EXCEPTABLE_ERROR);
-   EXPECT_NEAR(x, 0.0, EXCEPTABLE_ERROR);
-
-   integ.load();
-   integ.step();
-   integ.unload();
-
-   x = integ.getIndyVar();
-   EXPECT_NEAR(state[0],   (4.0/10.0), EXCEPTABLE_ERROR);
-   EXPECT_NEAR(state[1],  (-4.0/10.0), EXCEPTABLE_ERROR);
-   EXPECT_NEAR(state[2],  (M_PI/10.0), EXCEPTABLE_ERROR);
-   EXPECT_NEAR(state[3], (-M_PI/10.0), EXCEPTABLE_ERROR);
-   EXPECT_NEAR(x, 0.1, EXCEPTABLE_ERROR);
+    SA::Integrator  integ( h, udata);
+    integ.setIndyVar(4.0);
+
+    integ.integrate();
+    x = integ.getIndyVar();
+    EXPECT_NEAR(x, 4.1, EXCEPTABLE_ERROR);
+
+    integ.undo_integrate();
+    x = integ.getIndyVar();
+    EXPECT_NEAR(x, 4.0, EXCEPTABLE_ERROR);
 }