diff --git a/trick_source/trick_utils/SAIntegrator/README.md b/trick_source/trick_utils/SAIntegrator/README.md
index 3159a145..15e74572 100644
--- a/trick_source/trick_utils/SAIntegrator/README.md
+++ b/trick_source/trick_utils/SAIntegrator/README.md
@@ -13,6 +13,7 @@
 * [class RK2Integrator](#class-RK2Integrator)
 * [class RK4Integrator](#class-RK4Integrator)
 * [class RK3_8Integrator](#class-RK3_8Integrator)
+* [typedef Derivs2Func](#typedef-Derivs2Func)
 * [class EulerCromerIntegrator](#class-EulerCromerIntegrator)
 * [class ABM2Integrator](#class-ABM2Integrator)
 * [class ABM4Integrator](#class-ABM4Integrator)
@@ -53,10 +54,13 @@ This base-class represents a numerical **integrator**.
 |h        |```double```| Default integration step-size |
 |udata    |```void*``` | A pointer to user defined data that will be passed to user-defined functions when called by the Integrator. |
 
+
 ### Destructor
 
 #### ```virtual ~Integrator() {}```
 
+
+
 ### Public Member Functions
 
 <a id=method-Integrator::step></a>
@@ -172,6 +176,15 @@ where:
 | derivs_func |[```DerivsFunc```](#typedef-DerivsFunc)| Function thats generates the function (the derivatives) to be integrated. |
 |user_data    |```void*```  | A pointer to user defined data that will be passed to a DerivsFunc when called by the Integrator. |
 
+In addition to the above constructor, this class provides:
+
+* a copy constructor,
+* a destructor,
+* an assignment operator,
+* an insertion operator,
+* the public member functions inherited from [class Integrator](#class-Integrator),
+* and the following public member functions:
+
 ### Public Member Functions
 
 <a id=method-FirstOrderODEIntegrator::load></a>
@@ -307,6 +320,25 @@ Those inherited from [FirstOrderODEIntegrator](#class-FirstOrderODEIntegrator) p
 
 ### Constructor
 
+```
+FirstOrderODEVariableStepIntegrator( double h,
+                                     unsigned int N,
+                                     double* in_vars[],
+                                     double* out_vars[],
+                                     DerivsFunc dfunc,
+                                     void* udata);
+```
+[Constructor Parameters](#FOODEConstructorParameters) are those of [FirstOrderODEIntegrator](#class-FirstOrderODEIntegrator).
+
+In addition to the above constructor, this class provides:
+
+* a copy constructor,
+* a destructor,
+* an assignment operator,
+* an insertion operator,
+* the public member functions inherited from [FirstOrderODEIntegrator](#class-FirstOrderODEIntegrator),
+* and the following public member functions:
+
 <a id=FirstOrderODEVariableStepIntegrator::publicMemberFunctions></a>
 ### Public Member Functions
 
@@ -401,6 +433,15 @@ EulerIntegrator( double h,
 ```
 Constructor Parameters are those of [FirstOrderODEVariableStepIntegrator](#class-FirstOrderODEVariableStepIntegrator).
 
+In addition to the above constructor, this class provides:
+
+* a copy constructor,
+* a destructor,
+* an assignment operator,
+* an insertion operator,
+* the public member functions inherited from [FirstOrderODEVariableStepIntegrator](#class-FirstOrderODEVariableStepIntegrator),
+* and the following public member functions:
+
 ### Public Member Functions
 
 * All of the [Public Member Functions of FirstOrderODEVariableStepIntegrator](#FirstOrderODEVariableStepIntegrator::publicMemberFunctions), plus :
@@ -434,6 +475,15 @@ HeunsMethod( double h,
 ```
 [Constructor Parameters](#FOODEConstructorParameters) are those of [FirstOrderODEIntegrator](#class-FirstOrderODEIntegrator).
 
+In addition to the above constructor, this class provides:
+
+* a copy constructor,
+* a destructor,
+* an assignment operator,
+* an insertion operator,
+* the public member functions inherited from [FirstOrderODEVariableStepIntegrator](#class-FirstOrderODEVariableStepIntegrator),
+* and the following public member functions:
+
 ### Public Member Functions
 
 * All of the [Public Member Functions of FirstOrderODEVariableStepIntegrator](#FirstOrderODEVariableStepIntegrator::publicMemberFunctions).
@@ -463,10 +513,19 @@ RK2Integrator( double h,
                double* in_vars[],
                double* out_vars[],
                DerivsFunc func,
-           void* user_data)
+               void* user_data)
 ```
 [Constructor Parameters](#FOODEConstructorParameters) are those of [FirstOrderODEIntegrator](#class-FirstOrderODEIntegrator).
 
+In addition to the above constructor, this class provides:
+
+* a copy constructor,
+* a destructor,
+* an assignment operator,
+* an insertion operator,
+* the public member functions inherited from [FirstOrderODEVariableStepIntegrator](#class-FirstOrderODEVariableStepIntegrator),
+* and the following public member functions:
+
 ### Public Member Functions
 
 * All of the [Public Member Functions of FirstOrderODEVariableStepIntegrator](#FirstOrderODEVariableStepIntegrator::publicMemberFunctions).
@@ -501,6 +560,15 @@ RK4Integrator( double h,
 ```
 [Constructor Parameters](#FOODEConstructorParameters) are those of [FirstOrderODEIntegrator](#class-FirstOrderODEIntegrator).
 
+In addition to the above constructor, this class provides:
+
+* a copy constructor,
+* a destructor,
+* an assignment operator,
+* an insertion operator,
+* the public member functions inherited from [FirstOrderODEVariableStepIntegrator](#class-FirstOrderODEVariableStepIntegrator),
+* and the following public member functions:
+
 ### Public Member Functions
 
 * All of the [Public Member Functions of FirstOrderODEVariableStepIntegrator](#FirstOrderODEVariableStepIntegrator::publicMemberFunctions).
@@ -535,6 +603,15 @@ RK3_8Integrator( double h,
 [Constructor Parameters](#FOODEConstructorParameters) are those of
 [FirstOrderODEIntegrator](#class-FirstOrderODEIntegrator).
 
+In addition to the above constructor, this class provides:
+
+* a copy constructor,
+* a destructor,
+* an assignment operator,
+* an insertion operator,
+* the public member functions inherited from [FirstOrderODEVariableStepIntegrator](#class-FirstOrderODEVariableStepIntegrator),
+* and the following public member functions:
+
 ### Public Member Functions
 
 * All of the [Public Member Functions of FirstOrderODEVariableStepIntegrator](#FirstOrderODEVariableStepIntegrator::publicMemberFunctions).
@@ -567,6 +644,14 @@ ABM2Integrator ( double h,
 ```
 [Constructor Parameters](#FOODEConstructorParameters) are those of [FirstOrderODEIntegrator](#class-FirstOrderODEIntegrator).
 
+In addition to the above constructor, this class provides:
+
+* a copy constructor,
+* a destructor,
+* an assignment operator,
+* an insertion operator,
+* the public member functions inherited from [[FirstOrderODEIntegrator](#class-FirstOrderODEIntegrator).
+
 <a id=class-ABM4Integrator></a>
 ## class ABM4Integrator
 Derived from [FirstOrderODEIntegrator](#class-FirstOrderODEIntegrator).
@@ -589,6 +674,44 @@ ABM4Integrator ( double h,
 
 [Constructor Parameters](#FOODEConstructorParameters) are those of [FirstOrderODEIntegrator](#class-FirstOrderODEIntegrator).
 
+In addition to the above constructor, this class provides:
+
+* a copy constructor,
+* a destructor,
+* an assignment operator,
+* an insertion operator,
+* the public member functions inherited from [FirstOrderODEIntegrator](#class-FirstOrderODEIntegrator).
+
+
+<a id=typedef-Derivs2Func></a>
+## typedef Derivs2Func
+
+### Description
+This typedef defines a type of C/C++ function whose purpose is to populate
+an array of accelerations, given velocities and positions.
+
+```
+typedef void (*Derivs2Func)( double t, double x[], double v[], double a[], void* udata);
+```
+where:
+
+|Parameter|Type         |Direction|Description|
+|---------|-------------|---------|-----------|
+|t        |```double``` |IN       |Independent variable.|
+|x        |```double*```|IN       |Array of position values.|
+|v        |```double*```|IN       |Array of velocity values.|
+|a        |```double*```|OUT      |Array into which accelerations are to be returned.|
+|udata    |```void*```  |IN       |Pointer to user_data.|
+
+#### Example
+```
+void G( double t, double x[], double v[], double g_out[], void* udata) {
+    MassSpringDamper* msd = (MassSpringDamper*)udata;
+    g_out[0] = -(msd->k/msd->mass) * x[0]
+               -(msd->c/msd->mass) * v[0];
+}
+```
+
 <a id=class-EulerCromerIntegrator></a>
 ## class EulerCromerIntegrator
 Derived from [Integrator](#class-Integrator).
@@ -596,6 +719,12 @@ Derived from [Integrator](#class-Integrator).
 ### Description
 EulerCromer is integration method that conserves energy in oscillatory systems better than Runge-Kutta. So, it's good for mass-spring-damper systems, and orbital systems.
 
+It calculates the next state, from the current state as follows:
+
+![EulerCromerEqs](images/EulerCromerEqs.png)
+
+**a(v(n), x(n), t)** [above] is the function of type [```Derivs2Func```](#typedef-Derivs2Func) below.
+
 ### Data Members
 Those inherited from [Integrator](#class-Integrator) plus:
 
@@ -609,9 +738,9 @@ Those inherited from [Integrator](#class-Integrator) plus:
 | pos_out   |```double*```     |Protected|Position output array.|
 | vel_out   |```double*```     |Protected|Velocity output array.|
 | g_out     |```double*```     |Protected|Array of accelerations returned from gderivs.|
-| f_out     |```double*```     |Protected|Array of velocities returned from fderivs.|
-| gderivs   |[```DerivsFunc```](#typedef-DerivsFunc)|Protected|A function that returns accelerations.|
-| fderivs   |[```DerivsFunc```](#typedef-DerivsFunc)|Protected|A function that returns velocities.|
+| gderivs   |[```Derivs2Func```](#typedef-Derivs2Func)|Protected|A function that returns accelerations.|
+| last_h    |```double```|Value of h used in the last integration step.|
+
 
 ### Constructor
 ```
@@ -619,8 +748,7 @@ EulerCromerIntegrator(double dt,
                       int N,
                       double* xp[],
                       double* vp[],
-                      DerivsFunc gfunc,
-                      DerivsFunc ffunc,
+                      Derivs2Func gfunc,
                       void* user_data)
 ```
 
@@ -630,10 +758,17 @@ EulerCromerIntegrator(double dt,
 | N         |```int```    |Sets nDimensions above.|
 | xp        |```double*```|Sets pos_p above.|
 | vp        |```double*```|Sets vel_p above.|
-| gfunc |[```DerivsFunc```](#typedef-DerivsFunc)| Sets gderivs above. |
-| ffunc |[```DerivsFunc```](#typedef-DerivsFunc)| Sets fderivs above. |
+| gfunc |[```Derivs2Func```](#typedef-Derivs2Func)| Sets gderivs above. |
 |user_data  |```void*```  | Sets Integrator::user_data. |
 
+In addition to the above constructor, this class provides:
+
+* a copy constructor,
+* a destructor,
+* an assignment operator,
+* an insertion operator,
+* the public member functions inherited from [Integrator](#class-Integrator).
+
 ### Public Member Functions
 
 #### ```void step( double dt)```
diff --git a/trick_source/trick_utils/SAIntegrator/examples/CannonBall/CannonBall.cpp b/trick_source/trick_utils/SAIntegrator/examples/CannonBall/CannonBall.cpp
index 05a0d80a..a896b8a3 100644
--- a/trick_source/trick_utils/SAIntegrator/examples/CannonBall/CannonBall.cpp
+++ b/trick_source/trick_utils/SAIntegrator/examples/CannonBall/CannonBall.cpp
@@ -48,10 +48,7 @@ int main ( int argc, char* argv[]) {
     double dt = 0.01;
     double t = 0.0;
     SA::RK2Integrator integ(dt, 6, state_var_p, state_var_p, calc_derivs, NULL);
-
-    RootFinder root_finder(0.00000000001, Negative);
-    integ.add_Rootfinder(&root_finder, &impact);
-
+    integ.add_Rootfinder(0.00000000001, Negative, &impact);
     init_state(cannon);
     print_header();
     print_state( t, cannon);
diff --git a/trick_source/trick_utils/SAIntegrator/examples/MassSpringDamper/MassSpringDamper.cpp b/trick_source/trick_utils/SAIntegrator/examples/MassSpringDamper/MassSpringDamper.cpp
index a1aa2fc8..e207b6f0 100644
--- a/trick_source/trick_utils/SAIntegrator/examples/MassSpringDamper/MassSpringDamper.cpp
+++ b/trick_source/trick_utils/SAIntegrator/examples/MassSpringDamper/MassSpringDamper.cpp
@@ -13,7 +13,7 @@ void init_state ( MassSpringDamper& msd ) {
     msd.pos = 2.0;
     msd.vel = 0.0;
     msd.k = 1.0;
-    msd.c = 0.0;
+    msd.c = 0.2;
     msd.mass = 5.0;
 }
 void print_header() {
@@ -22,31 +22,24 @@ void print_header() {
 void print_state( double t, MassSpringDamper& msd ) {
     printf ("%10.10f, %10.10f, %10.10f\n", t, msd.pos, msd.vel);
 }
-void G( double t, double x[], double g_out[], void* udata) {
+void G( double t, double x[], double v[], double g_out[], void* udata) {
     MassSpringDamper* msd = (MassSpringDamper*)udata;
-    g_out[0] = -(msd->k/msd->mass) * x[0];
-}
-void F( double t, double v[], double f_out[], void* udata) {
-    f_out[0] = v[0];
-    f_out[1] = v[1];
+    g_out[0] = -(msd->k/msd->mass) * x[0]
+               -(msd->c/msd->mass) * v[0];
 }
 int main ( int argc, char* argv[]) {
     MassSpringDamper msd;
+    double time = 0.0;
     double* x_p[N_DIMENSIONS] = { &msd.pos };
     double* v_p[N_DIMENSIONS] = { &msd.vel };
-    unsigned count = 0;
     double dt = 0.001;
-    double time = count * dt;
     init_state(msd);
     print_header();
     print_state(time, msd);
-    SA::EulerCromerIntegrator integ(dt, N_DIMENSIONS, x_p, v_p, G, F, &msd);
+    SA::EulerCromerIntegrator integ(dt, N_DIMENSIONS, x_p, v_p, G, &msd);
     while (time < 50) {
-        integ.load();
-        integ.step();
-        integ.unload();
-        count ++;
-        time = count * dt;
+        integ.integrate();
+        time = integ.getIndyVar();
         print_state(time, msd);
     }
 }
diff --git a/trick_source/trick_utils/SAIntegrator/examples/MassSpringDamper/plot_position.py b/trick_source/trick_utils/SAIntegrator/examples/MassSpringDamper/plot_position.py
old mode 100644
new mode 100755
diff --git a/trick_source/trick_utils/SAIntegrator/examples/Orbit/Orbit.cpp b/trick_source/trick_utils/SAIntegrator/examples/Orbit/Orbit.cpp
index 85c56f0b..e68c7d29 100644
--- a/trick_source/trick_utils/SAIntegrator/examples/Orbit/Orbit.cpp
+++ b/trick_source/trick_utils/SAIntegrator/examples/Orbit/Orbit.cpp
@@ -25,33 +25,25 @@ void print_state( double t, Orbit& orbit ) {
     printf ("%10.10f, %10.10f, %10.10f, %10.10f, %10.10f\n",
              t, orbit.pos[0], orbit.pos[1], orbit.vel[0], orbit.vel[1]);
 }
-void G( double t, double x[], double g_out[], void* udata) {
+void G( double t, double pos[], double vel[], double acc_out[], void* udata) {
     Orbit* orbit = (Orbit*)udata;
-    double d = sqrt( x[0]*x[0] + x[1]*x[1]);
-    g_out[0] = -x[0] * GRAVITATIONAL_CONSTANT * orbit->planet_mass / (d*d*d);
-    g_out[1] = -x[1] * GRAVITATIONAL_CONSTANT * orbit->planet_mass / (d*d*d);
-}
-void F( double t, double v[], double f_out[], void* udata) {
-    f_out[0] = v[0];
-    f_out[1] = v[1];
+    double d = sqrt( pos[0]*pos[0] + pos[1]*pos[1]);
+    acc_out[0] = -pos[0] * GRAVITATIONAL_CONSTANT * orbit->planet_mass / (d*d*d);
+    acc_out[1] = -pos[1] * GRAVITATIONAL_CONSTANT * orbit->planet_mass / (d*d*d);
 }
 int main ( int argc, char* argv[]) {
     Orbit orbit;
+    double time = 0.0;
     double* x_p[2] = { &orbit.pos[0], &orbit.pos[1] };
     double* v_p[2] = { &orbit.vel[0], &orbit.vel[1] };
-    unsigned count = 0;
     double dt = 0.1;
-    double time = count * dt;
     init_state(orbit);
     print_header();
     print_state(time, orbit);
-    SA::EulerCromerIntegrator integ(dt, 2, x_p, v_p, G, F, &orbit);
-    while (time < 6000) {
-        integ.load();
-        integ.step();
-        integ.unload();
-        count ++;
-        time = count * dt;
+    SA::EulerCromerIntegrator integ(dt, 2, x_p, v_p, G, &orbit);
+    while (time < 5600) {
+        integ.integrate();
+        time = integ.getIndyVar();
         print_state(time, orbit);
     }
 }
diff --git a/trick_source/trick_utils/SAIntegrator/images/EulerCromerEqs.png b/trick_source/trick_utils/SAIntegrator/images/EulerCromerEqs.png
new file mode 100644
index 00000000..dfa52e31
Binary files /dev/null and b/trick_source/trick_utils/SAIntegrator/images/EulerCromerEqs.png differ
diff --git a/trick_source/trick_utils/SAIntegrator/include/Rootfinder.hh b/trick_source/trick_utils/SAIntegrator/include/Rootfinder.hh
index 43e01425..327d0e0c 100644
--- a/trick_source/trick_utils/SAIntegrator/include/Rootfinder.hh
+++ b/trick_source/trick_utils/SAIntegrator/include/Rootfinder.hh
@@ -1,3 +1,4 @@
+#include <iostream>
 
 typedef enum {
     Negative = -1,
@@ -12,6 +13,7 @@ class RootFinder {
         RootFinder();
         RootFinder (double tolerance, SlopeConstraint constraint);
         double find_roots( double x, double f_error );
+        friend std::ostream& operator<<(std::ostream& os, const RootFinder& rf);
     private:
         double f_upper;
         double x_upper;
diff --git a/trick_source/trick_utils/SAIntegrator/include/SAIntegrator.hh b/trick_source/trick_utils/SAIntegrator/include/SAIntegrator.hh
index 30153ab1..41e148c4 100644
--- a/trick_source/trick_utils/SAIntegrator/include/SAIntegrator.hh
+++ b/trick_source/trick_utils/SAIntegrator/include/SAIntegrator.hh
@@ -1,4 +1,6 @@
 #include "Rootfinder.hh"
+#include <iostream>
+
 namespace SA {
 
     class Integrator {
@@ -9,7 +11,8 @@ namespace SA {
         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) {};
+        Integrator();
+        Integrator(double h, void* udata);
         virtual ~Integrator() {}
         virtual void load();
         virtual void unload();
@@ -18,8 +21,11 @@ namespace SA {
         virtual double undo_integrate();
         double getIndyVar();
         void setIndyVar(double v);
+        friend std::ostream& operator<<(std::ostream& os, const SA::Integrator& I);
     };
 
+    std::ostream& operator<<(std::ostream& os, const Integrator& I);
+
     typedef void (*DerivsFunc)( double x, double state[], double derivs[], void* udata);
 
     class FirstOrderODEIntegrator : public Integrator {
@@ -30,21 +36,23 @@ 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.
-     // void advanceIndyVar();   // Inherited from Integrator::advanceIndyVar()
     public:
-        FirstOrderODEIntegrator( double h, int N, double* in_vars[], double* out_vars[], DerivsFunc dfunc, void* udata);
+        FirstOrderODEIntegrator();
+        FirstOrderODEIntegrator( double h, unsigned int N, double* in_vars[], double* out_vars[], DerivsFunc dfunc, void* udata);
+        FirstOrderODEIntegrator(const FirstOrderODEIntegrator& other);
+        FirstOrderODEIntegrator& operator=( const FirstOrderODEIntegrator& other);
         virtual ~FirstOrderODEIntegrator();
         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);
+
+        friend std::ostream& operator<<(std::ostream& os, const SA::FirstOrderODEIntegrator& I);
     };
+    std::ostream& operator<<(std::ostream& os, const FirstOrderODEIntegrator& I);
 
     typedef double (*RootErrorFunc)( double x, double state[], RootFinder* root_finder, void* udata);
 
@@ -54,82 +62,96 @@ namespace SA {
         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();
+        FirstOrderODEVariableStepIntegrator( double h, unsigned int N, double* in_vars[], double* out_vars[], DerivsFunc dfunc, void* udata);
+        FirstOrderODEVariableStepIntegrator( const FirstOrderODEVariableStepIntegrator& other);
+        FirstOrderODEVariableStepIntegrator& operator=( const FirstOrderODEVariableStepIntegrator& rhs);
         ~FirstOrderODEVariableStepIntegrator();
-        // 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
+        void add_Rootfinder( double tolerance, SlopeConstraint constraint, RootErrorFunc rfunc);
+        friend std::ostream& operator<<(std::ostream& os, const SA::FirstOrderODEVariableStepIntegrator& I);
     private:
         void find_roots(double h, unsigned int depth);
     };
+    std::ostream& operator<<(std::ostream& os, const FirstOrderODEVariableStepIntegrator& I);
 
     class EulerIntegrator : public FirstOrderODEVariableStepIntegrator {
         public:
         double *derivs;
-        EulerIntegrator( double h, int N, double* in_vars[], double* out_vars[], DerivsFunc dfunc, void* udata);
+        EulerIntegrator();
+        EulerIntegrator( double h, unsigned int N, double* in_vars[], double* out_vars[], DerivsFunc dfunc, void* udata);
+        EulerIntegrator(const EulerIntegrator& other);
+        EulerIntegrator& operator=( const EulerIntegrator& rhs);
         ~EulerIntegrator();
-        // 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
+        void variable_step( double h);
+        friend std::ostream& operator<<(std::ostream& os, const SA::EulerIntegrator& I);
     };
+    std::ostream& operator<<(std::ostream& os, const EulerIntegrator& I);
 
     class HeunsMethod : public FirstOrderODEVariableStepIntegrator {
         public:
         double *wstate;
         double *derivs[2];
-        HeunsMethod( double h, int N, double* in_vars[], double* out_vars[], DerivsFunc dfunc, void* udata);
+        HeunsMethod();
+        HeunsMethod( double h, unsigned int N, double* in_vars[], double* out_vars[], DerivsFunc dfunc, void* udata);
+        HeunsMethod(const HeunsMethod& other);
+        HeunsMethod& operator=( const HeunsMethod& rhs);
         ~HeunsMethod();
-        void variable_step( double h); // Overrides FirstOrderODEVariableStepIntegrator::variable_step()
+        void variable_step( double h);
+        friend std::ostream& operator<<(std::ostream& os, const SA::HeunsMethod& I);
     };
+    std::ostream& operator<<(std::ostream& os, const HeunsMethod& I);
 
     class RK2Integrator : public FirstOrderODEVariableStepIntegrator {
         public:
         double *wstate;
         double *derivs[2];
-        RK2Integrator( double h, int N, double* in_vars[], double* out_vars[], DerivsFunc derivs_func, void* udata);
+        RK2Integrator();
+        RK2Integrator( double h, unsigned int N, double* in_vars[], double* out_vars[], DerivsFunc derivs_func, void* udata);
+        RK2Integrator(const RK2Integrator& other);
+        RK2Integrator& operator=( const RK2Integrator& rhs);
         ~RK2Integrator();
-        void variable_step( double h); // Overrides FirstOrderODEVariableStepIntegrator::variable_step()
+        void variable_step( double h);
+        friend std::ostream& operator<<(std::ostream& os, const SA::RK2Integrator& I);
     };
+    std::ostream& operator<<(std::ostream& os, const RK2Integrator& I);
 
     class RK4Integrator : public FirstOrderODEVariableStepIntegrator {
         public:
         double *wstate[3];
         double *derivs[4];
-        RK4Integrator( double h, int N, double* in_vars[], double* out_vars[], DerivsFunc dfunc, void* udata);
+
+        RK4Integrator();
+        RK4Integrator( double h, unsigned int N, double* in_vars[], double* out_vars[], DerivsFunc dfunc, void* udata);
+        RK4Integrator(const RK4Integrator& other);
+        RK4Integrator& operator=( const RK4Integrator& rhs);
         ~RK4Integrator();
-        void variable_step( double h); // Overrides FirstOrderODEVariableStepIntegrator::variable_step()
+        void variable_step( double h);
+        friend std::ostream& operator<<(std::ostream& os, const SA::RK4Integrator& I);
     };
+    std::ostream& operator<<(std::ostream& os, const RK4Integrator& I);
 
     class RK3_8Integrator : public FirstOrderODEVariableStepIntegrator {
         public:
-        double *wstate[4];
+        double *wstate[3];
         double *derivs[4];
-        RK3_8Integrator( double h, int N, double* in_vars[], double* out_vars[], DerivsFunc dfunc, void* udata);
+
+        RK3_8Integrator();
+        RK3_8Integrator( double h, unsigned int N, double* in_vars[], double* out_vars[], DerivsFunc dfunc, void* udata);
+        RK3_8Integrator(const RK3_8Integrator& other);
+        RK3_8Integrator& operator=( const RK3_8Integrator& rhs);
         ~RK3_8Integrator();
-        void variable_step( double h); // Overrides FirstOrderODEVariableStepIntegrator::variable_step()
+        void variable_step( double h);
+
+        friend std::ostream& operator<<(std::ostream& os, const SA::RK3_8Integrator& I);
     };
+    std::ostream& operator<<(std::ostream& os, const RK3_8Integrator& I);
+
+    typedef void (*Derivs2Func)( double t, double x[], double v[], double derivs[], void* udata);
 
     class EulerCromerIntegrator : public Integrator {
     protected:
@@ -141,19 +163,24 @@ namespace SA {
         double* pos_out;
         double* vel_out;
         double* g_out;
-        double* f_out;
-        DerivsFunc gderivs;
-        DerivsFunc fderivs;
+        double last_h;
+        Derivs2Func gderivs;
 
     public:
-        EulerCromerIntegrator(double dt, int N, double* xp[], double* vp[], DerivsFunc g, DerivsFunc f , void* udata);
+        EulerCromerIntegrator();
+        EulerCromerIntegrator(double dt, unsigned int N, double* xp[], double* vp[], Derivs2Func g, void* udata);
+        EulerCromerIntegrator(const EulerCromerIntegrator& other);
+        EulerCromerIntegrator& operator=( const EulerCromerIntegrator& rhs);
         ~EulerCromerIntegrator();
+        void advanceIndyVar(double h);
         void step( double dt);
         void step();
         void load();
         void unload();
         double undo_integrate();
+        friend std::ostream& operator<<(std::ostream& os, const SA::EulerCromerIntegrator& I);
     };
+    std::ostream& operator<<(std::ostream& os, const EulerCromerIntegrator& I);
 
     // AdamsBashforthMoulton Self Priming
     class ABMIntegrator : public FirstOrderODEIntegrator {
@@ -165,11 +192,15 @@ namespace SA {
         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(unsigned int history_size, double h, unsigned int N, double* in_vars[], double* out_vars[], DerivsFunc dfunc, void* udata);
+        ABMIntegrator(const ABMIntegrator& other);
+        ABMIntegrator& operator=( const ABMIntegrator& rhs);
         ~ABMIntegrator();
         void step();
         double undo_integrate();
+        friend std::ostream& operator<<(std::ostream& os, const ABMIntegrator& I);
     };
+    std::ostream& operator<<(std::ostream& os, const ABMIntegrator& I);
 
     // AdamsBashforthMoulton 2 with RK2 Priming
     class ABM2Integrator : public FirstOrderODEIntegrator {
@@ -179,13 +210,17 @@ namespace SA {
         unsigned int bufferSize;
         unsigned int hix;
         unsigned int currentHistorySize;
-        FirstOrderODEVariableStepIntegrator* priming_integrator;
+        RK2Integrator* priming_integrator;
     public:
-        ABM2Integrator( double h, int N, double* in_vars[], double* out_vars[], DerivsFunc dfunc, void* udata);
+        ABM2Integrator( double h, unsigned int N, double* in_vars[], double* out_vars[], DerivsFunc dfunc, void* udata);
+        ABM2Integrator(const ABM2Integrator& other);
+        ABM2Integrator& operator=( const ABM2Integrator& rhs);
         ~ABM2Integrator();
         void step();
         double undo_integrate();
+        friend std::ostream& operator<<(std::ostream& os, const ABM2Integrator& I);
     };
+    std::ostream& operator<<(std::ostream& os, const ABM2Integrator& I);
 
     // AdamsBashforthMoulton 4 with RK4 Priming
     class ABM4Integrator : public FirstOrderODEIntegrator {
@@ -195,11 +230,15 @@ namespace SA {
         unsigned int bufferSize;
         unsigned int hix;
         unsigned int currentHistorySize;
-        FirstOrderODEVariableStepIntegrator* priming_integrator;
+        RK4Integrator* priming_integrator;
     public:
-        ABM4Integrator( double h, int N, double* in_vars[], double* out_vars[], DerivsFunc dfunc, void* udata);
+        ABM4Integrator( double h, unsigned int N, double* in_vars[], double* out_vars[], DerivsFunc dfunc, void* udata);
+        ABM4Integrator(const ABM4Integrator& other);
+        ABM4Integrator& operator=( const ABM4Integrator& rhs);
         ~ABM4Integrator();
         void step();
         double undo_integrate();
+        friend std::ostream& operator<<(std::ostream& os, const ABM4Integrator& I);
     };
+    std::ostream& operator<<(std::ostream& os, const ABM4Integrator& I);
 }
diff --git a/trick_source/trick_utils/SAIntegrator/src/RootFinder.cpp b/trick_source/trick_utils/SAIntegrator/src/RootFinder.cpp
index d3044dee..3fd0f62b 100644
--- a/trick_source/trick_utils/SAIntegrator/src/RootFinder.cpp
+++ b/trick_source/trick_utils/SAIntegrator/src/RootFinder.cpp
@@ -14,12 +14,18 @@ void RootFinder::init( double tolerance, SlopeConstraint constraint) {
 void RootFinder::init() {
     init(0.00000000001, Unconstrained);
 }
-RootFinder::RootFinder (double tolerance, SlopeConstraint constraint) {
-    init(tolerance, constraint);
-}
+// Default Constructor
 RootFinder::RootFinder () {
     init();
 }
+// Constructor
+RootFinder::RootFinder (double tolerance, SlopeConstraint constraint) {
+    init(tolerance, constraint);
+}
+// Copy Constructor: The default copy_constructor should suffice.
+// Assignment Operator: The default assignment operator should suffice.
+// Destructor: The default destructor should suffice.
+
 // Given the values of the independent variable, and the function,
 // estimate the distance of the independent variable from functions root.
 double RootFinder::find_roots( double x, double f_error ) {
@@ -86,3 +92,15 @@ double RootFinder::find_roots( double x, double f_error ) {
     iterations = 0;
     return (DBL_MAX);
 }
+
+std::ostream& operator<<(std::ostream& os, const RootFinder& rf) {
+    os << "\n--- RootFinder ---";
+    os << "\nupper: {" << rf.f_upper << ',' << rf.x_upper << ',' << rf.upper_set << "}";
+    os << "\nlower: {" << rf.f_lower << ',' << rf.x_lower << ',' << rf.lower_set << "}";
+    os << "\nprev_f_error: " << rf.prev_f_error;
+    os << "\nf_error_tol:" << rf.prev_f_error;
+    os << "\niterations:" << rf.prev_f_error;
+    os << "\nslope_constraint:" << rf.slope_constraint;
+    os << "\nf_slope:" << rf.f_slope;
+    return os;
+}
diff --git a/trick_source/trick_utils/SAIntegrator/src/SAIntegrator.cpp b/trick_source/trick_utils/SAIntegrator/src/SAIntegrator.cpp
index f93ecedb..3232c65d 100644
--- a/trick_source/trick_utils/SAIntegrator/src/SAIntegrator.cpp
+++ b/trick_source/trick_utils/SAIntegrator/src/SAIntegrator.cpp
@@ -1,8 +1,33 @@
 
 #include <stdlib.h>
-#include <iostream>
+#include <stdexcept>
+#include <iostream>  // std::cout, std::cerr
+#include <algorithm> // std::copy
 #include "SAIntegrator.hh"
 
+static void nil_dfunc( double x,
+                       double state[] __attribute__((unused)),
+                       double derivs[],
+                       void*  udata __attribute__((unused))) {
+    derivs[0] = 0.0;
+}
+static void nil_gfunc( double t,
+                       double x[] __attribute__((unused)),
+                       double v[] __attribute__((unused)),
+                       double derivs[],
+                       void*  udata __attribute__((unused))) {
+    derivs[0] = 0.0;
+}
+// ------------------------------------------------------------
+// Class Integrator
+// ------------------------------------------------------------
+// Default Constructor
+SA::Integrator::Integrator()
+: X_in(0.0), X_out(0.0), default_h(1.0), user_data(NULL) {};
+// Constructor
+SA::Integrator::Integrator(double h, void* udata)
+: X_in(0.0), X_out(0.0), default_h(h), user_data(udata) {};
+
 void SA::Integrator::advanceIndyVar() {
     X_out = X_in + default_h;
 }
@@ -30,21 +55,80 @@ double SA::Integrator::getIndyVar() {
 void SA::Integrator::setIndyVar(double v) {
     X_out = v;
 }
+// Insertion Operator
+std::ostream& SA::operator<<(std::ostream& os, const Integrator& I) {
+    os << "\n--- Integrator ---";
+    os << "\nX_in      : " << I.X_in;
+    os << "\nX_out     : " << I.X_out;
+    os << "\ndefault_h : " << I.default_h;
+    os << "\nuser_data : " << I.user_data;
+    return os;
+}
 
 // ------------------------------------------------------------
+// Class FirstOrderODEIntegrator
+// ------------------------------------------------------------
+// Default Constructor
+SA::FirstOrderODEIntegrator::FirstOrderODEIntegrator()
+: Integrator() {
+    state_size = 0;
+    inVars = NULL;
+    outVars = NULL;
+    inState = new double[state_size];
+    outState = new double[state_size];
+    derivs_func = nil_dfunc;
+}
+// Constructor
 SA::FirstOrderODEIntegrator::FirstOrderODEIntegrator(
-    double h, int N, double* in_vars[], double* out_vars[], DerivsFunc dfunc, void* udata)
+    double h, unsigned int N, double* in_vars[], double* out_vars[], DerivsFunc dfunc, void* udata)
 : Integrator(h, udata) {
     state_size = N;
     inVars = in_vars;
     outVars = out_vars;
     inState = new double[state_size];
     outState = new double[state_size];
+    if (dfunc == NULL) throw std::invalid_argument("dfunc must be non-NULL.");
     derivs_func = dfunc;
 }
+// Copy Constructor
+SA::FirstOrderODEIntegrator::FirstOrderODEIntegrator(const FirstOrderODEIntegrator& other)
+: Integrator(other) {
+    state_size = other.state_size;
+    inVars = other.inVars;
+    outVars = other.outVars;
+    inState = new double[state_size];
+    std::copy( other.inState, other.inState + state_size, inState);
+    outState = new double[state_size];
+    std::copy( other.outState, other.outState + state_size, outState);
+    derivs_func = other.derivs_func;
+}
+// Assignment Operator
+SA::FirstOrderODEIntegrator& SA::FirstOrderODEIntegrator::operator=( const SA::FirstOrderODEIntegrator& rhs) {
+    if (this != &rhs) {
+        // Call base assignment operator
+        Integrator::operator=(rhs);
+        // Duplicate rhs arrays
+        double* new_inState = new double[rhs.state_size];
+        std::copy( rhs.inState, rhs.inState + rhs.state_size, new_inState);
+        double* new_outState = new double[rhs.state_size];
+        std::copy( rhs.outState, rhs.outState + rhs.state_size, new_outState);
+        // Delete lhs arrays & replace with rhs arrays
+        if (inState != NULL) delete inState;
+        inState = new_inState;
+        if (outState != NULL) delete outState;
+        outState = new_outState;
+        // Copy primitive members
+        state_size = rhs.state_size;
+        inVars = rhs.inVars;
+        outVars = rhs.outVars;
+        derivs_func = rhs.derivs_func;
+    }
+    return *this;
+}
+// Destructor
 SA::FirstOrderODEIntegrator::~FirstOrderODEIntegrator() {
-    delete inState;
-    delete outState;
+    if (inState != NULL) delete inState;
+    if (outState != NULL) delete outState;
 }
 void SA::FirstOrderODEIntegrator::load() {
     if (inVars != NULL) {
@@ -82,13 +166,16 @@ double** SA::FirstOrderODEIntegrator::set_out_vars( double* out_vars[]){
     return (ret);
 }
 double SA::FirstOrderODEIntegrator::undo_integrate() {
-    for (unsigned int i=0 ; i<state_size; i++ ) {
-        *(inVars[i]) = inState[i];
+    if (inVars != NULL) {
+        for (unsigned int i=0 ; i<state_size; i++ ) {
+            *(inVars[i]) = inState[i];
+        }
+        std::copy(inState, inState + state_size, outState);
+        return (SA::Integrator::undo_integrate());
+    } else {
+        std::cerr << "Error: SA::FirstOrderODEIntegrator::load(). inVars is not set." << std::endl;
     }
-    for (unsigned int i=0 ; i<state_size; i++ ) {
-        outState[i] = inState[i];
-    }
-    return (SA::Integrator::undo_integrate());
+    return (0.0);
 }
 void SA::FirstOrderODEIntegrator::step() {
     double derivs[state_size];
@@ -98,24 +185,88 @@ void SA::FirstOrderODEIntegrator::step() {
     }
     advanceIndyVar();
 }
+static std::ostream& stream_double_array(std::ostream& os, unsigned int n, double* d_array) {
+    if (d_array == NULL) {
+        os << " NULL\n";
+    } else {
+        os << "[";
+        for (int i=0; i<n ; i++) {
+            if (i==0) { os << d_array[i]; }
+            else { os << " ," << d_array[i]; }
+        }
+        os << "]";
+    }
+    return os;
+}
+// Insertion Operator
+std::ostream& SA::operator<<(std::ostream& os, const FirstOrderODEIntegrator& I) {
+        os << (SA::Integrator)I;
+        os << "\n--- FirstOrderODEIntegrator ---";
+        os << "\nstate_size: " << I.state_size;
+        os << "\ninState   :";
+        stream_double_array(os, I.state_size, I.inState);
+        os << "\noutState  :";
+        stream_double_array(os, I.state_size, I.outState);
+        os << "\ninVars    : " << (void*)I.inVars;
+        os << "\noutVars   : " << (void*)I.outVars;
+        return os;
+}
 
 // ------------------------------------------------------------
+// Class FirstOrderODEVariableStepIntegrator
+// ------------------------------------------------------------
+// Default Constructor
+SA::FirstOrderODEVariableStepIntegrator::FirstOrderODEVariableStepIntegrator()
+: FirstOrderODEIntegrator() {
+    root_finder = (RootFinder*) NULL;
+    root_error_func = (RootErrorFunc) NULL;
+}
+// Constructor
 SA::FirstOrderODEVariableStepIntegrator::FirstOrderODEVariableStepIntegrator(
-    double h, int N, double* in_vars[], double* out_vars[], DerivsFunc dfunc, void* udata)
+    double h, unsigned int N, double* in_vars[], double* out_vars[], DerivsFunc dfunc, void* udata)
 : FirstOrderODEIntegrator(h, N, in_vars, out_vars, dfunc, udata) {
     root_finder = (RootFinder*) NULL;
-    root_error_func = (RootErrorFunc)NULL;
+    root_error_func = (RootErrorFunc) NULL;
+}
+// Copy Constructor
+SA::FirstOrderODEVariableStepIntegrator::FirstOrderODEVariableStepIntegrator(
+const FirstOrderODEVariableStepIntegrator& other)
+: FirstOrderODEIntegrator(other) {
+    if (other.root_finder != NULL) {
+        root_finder = new RootFinder(*(other.root_finder));
+    } else {
+        root_finder = NULL;
+    }
+    root_error_func = other.root_error_func;
+}
+// Assignment Operator
+SA::FirstOrderODEVariableStepIntegrator& SA::
+FirstOrderODEVariableStepIntegrator::operator=( const SA::FirstOrderODEVariableStepIntegrator& rhs) {
+    FirstOrderODEIntegrator::operator=(rhs);
+    if (root_finder != NULL) delete root_finder;
+    if (rhs.root_finder != NULL) {
+        root_finder = new RootFinder(*(rhs.root_finder));
+    } else {
+        root_finder = NULL;
+    }
+    root_error_func = rhs.root_error_func;
+    return *this;
+}
+// Destructor
+SA::FirstOrderODEVariableStepIntegrator::~FirstOrderODEVariableStepIntegrator() {
+    if (root_finder != NULL) {
+        delete root_finder;
+    }
 }
 double SA::FirstOrderODEVariableStepIntegrator::undo_integrate() {
     SA::FirstOrderODEIntegrator::undo_integrate();
     return (last_h);
 }
-SA::FirstOrderODEVariableStepIntegrator::~FirstOrderODEVariableStepIntegrator() {
-    // User is responsible for deleting the root_finder.
-}
-void SA::FirstOrderODEVariableStepIntegrator::add_Rootfinder(
-    RootFinder* rootFinder, RootErrorFunc rfunc) {
-    root_finder = rootFinder;
+void SA::FirstOrderODEVariableStepIntegrator::add_Rootfinder( double tolerance, SlopeConstraint constraint, RootErrorFunc rfunc) {
+    root_finder = new RootFinder(tolerance, constraint);
+    if (rfunc == NULL) {
+        throw std::invalid_argument("OOPS! RootErrorFunc function-pointer must be NULL.");
+    }
     root_error_func = rfunc;
 }
 void SA::FirstOrderODEVariableStepIntegrator::advanceIndyVar(double h) {
@@ -151,15 +302,63 @@ void SA::FirstOrderODEVariableStepIntegrator::step() {
         find_roots( default_h, 5 );
     }
 }
+// Insertion Operator
+std::ostream& SA::operator<<(std::ostream& os, const FirstOrderODEVariableStepIntegrator& I) {
+    os << (SA::FirstOrderODEIntegrator)I;
+    os << "\n--- FirstOrderODEVariableStepIntegrator ---";
+    if (I.root_error_func == NULL) {
+        os << "\nroot_error_func: NULL.";
+    } else {
+        os << "\nroot_error_func: " << (void*)I.root_error_func;
+    }
+    if (I.root_finder == NULL) {
+        os << "\nroot_finder: NULL.";
+    } else {
+        os << "\nroot_finder -> " << *(I.root_finder);
+    }
+    return os;
+}
 
 // ------------------------------------------------------------
+// Class EulerIntegrator
+// ------------------------------------------------------------
+// Default Constructor
+SA::EulerIntegrator::EulerIntegrator()
+: FirstOrderODEVariableStepIntegrator() {
+    derivs = new double[state_size];
+}
+// Constructor
 SA::EulerIntegrator::EulerIntegrator(
-    double h, int N, double* in_vars[], double* out_vars[], DerivsFunc dfunc, void* udata)
+    double h, unsigned int N, double* in_vars[], double* out_vars[], DerivsFunc dfunc, void* udata)
     : FirstOrderODEVariableStepIntegrator(h, N, in_vars, out_vars, dfunc, udata) {
     derivs = new double[N];
 }
+// Copy Constructor
+SA::EulerIntegrator::EulerIntegrator(const EulerIntegrator& other)
+    : FirstOrderODEVariableStepIntegrator(other) {
+    derivs = new double[state_size];
+    std::copy( other.derivs, other.derivs + state_size, derivs);
+}
+// Assignment Operator
+SA::EulerIntegrator& SA::EulerIntegrator::operator=( const SA::EulerIntegrator& rhs) {
+    if (this != &rhs) {
+        double* new_derivs = 0;
+        try {
+            new_derivs = new double[rhs.state_size];
+        } catch (...) {
+            delete new_derivs;
+            throw;
+        }
+        FirstOrderODEVariableStepIntegrator::operator=(rhs);
+        if (derivs != NULL) delete derivs;
+        derivs = new_derivs;
+        std::copy( rhs.derivs, rhs.derivs + state_size, derivs);
+    }
+    return *this;
+}
+//Destructor
 SA::EulerIntegrator::~EulerIntegrator() {
-    delete derivs;
+    if (derivs != NULL) delete derivs;
 }
 void SA::EulerIntegrator::variable_step( double h) {
     (*derivs_func)( X_in, inState, derivs, user_data);
@@ -168,19 +367,72 @@ void SA::EulerIntegrator::variable_step( double h) {
     }
     advanceIndyVar(h);
 }
+// Insertion Operator
+std::ostream& SA::operator<<(std::ostream& os, const EulerIntegrator& I) {
+    os << (SA::FirstOrderODEVariableStepIntegrator)I;
+    os << "\n--- EulerIntegrator ---";
+    os << "\nderivs    :";
+    stream_double_array(os, I.state_size, I.derivs);
+    return os;
+}
 // ------------------------------------------------------------
+// Class HeunsMethod
+// ------------------------------------------------------------
+// Default Constructor
+SA::HeunsMethod::HeunsMethod()
+: FirstOrderODEVariableStepIntegrator() {
+    wstate = new double[state_size];
+    for (unsigned int i = 0; i < 2 ; i++) {
+        derivs[i] = new double[state_size];
+    }
+}
+// Constructor
 SA::HeunsMethod::HeunsMethod(
-    double h, int N, double* in_vars[], double* out_vars[], DerivsFunc dfunc, void* udata)
+    double h, unsigned int N, double* in_vars[], double* out_vars[], DerivsFunc dfunc, void* udata)
     : FirstOrderODEVariableStepIntegrator(h, N, in_vars, out_vars, dfunc, udata) {
     wstate = new double[N];
     for (unsigned int i = 0; i < 2 ; i++) {
         derivs[i] = new double[N];
     }
 }
+// Copy Constructor
+SA::HeunsMethod::HeunsMethod(const HeunsMethod& other)
+    : FirstOrderODEVariableStepIntegrator(other) {
+    wstate = new double[state_size];
+    std::copy(other.wstate, other.wstate+state_size, wstate);
+    for (unsigned int i=0 ; i<2; i++ ) {
+        derivs[i] = new double[state_size];
+        std::copy(other.derivs[i], other.derivs[i]+state_size, derivs[i]);
+    }
+}
+// Assignment Operator
+SA::HeunsMethod& SA::HeunsMethod::operator=( const SA::HeunsMethod& rhs) {
+    if (this != &rhs) {
+        // Call base assignment operator
+        FirstOrderODEVariableStepIntegrator::operator=(rhs);
+        // Duplicate rhs arrays
+        double* new_wstate = new double[rhs.state_size];
+        std::copy(rhs.wstate, rhs.wstate + rhs.state_size, new_wstate);
+        double* new_derivs[2];
+        for (unsigned int i=0 ; i<2; i++ ) {
+            new_derivs[i] = new double[rhs.state_size];
+            std::copy(rhs.derivs[i], rhs.derivs[i] + rhs.state_size, new_derivs[i]);
+        }
+        // Delete lhs arrays & replace with rhs arrays
+        if (wstate != NULL) delete wstate;
+        wstate = new_wstate;
+        for (unsigned int i=0 ; i<2; i++ ) {
+            if (derivs[i] != NULL) delete derivs[i];
+            derivs[i] = new_derivs[i];
+        }
+    }
+    return *this;
+}
+// Destructor
 SA::HeunsMethod::~HeunsMethod() {
-    delete wstate;
+    if (wstate != NULL) delete wstate;
     for (unsigned int i = 0; i < 2 ; i++) {
-        delete derivs[i];
+        if (derivs[i] != NULL) delete derivs[i];
     }
 }
 void SA::HeunsMethod::variable_step( double h) {
@@ -194,19 +446,78 @@ void SA::HeunsMethod::variable_step( double h) {
     }
     advanceIndyVar(h);
 }
+// Insertion Operator
+std::ostream& SA::operator<<(std::ostream& os, const HeunsMethod& I) {
+    os << (SA::FirstOrderODEVariableStepIntegrator)I;
+    os << "\n--- HeunsMethod ---";
+    os << "\nderivs    :[";
+    for (int i=0; i<2 ; i++) {
+        os << "\n";
+        stream_double_array(os, I.state_size, I.derivs[i]);
+    }
+    os << "\n]";
+    os << "\nwstate    :";
+    stream_double_array(os, I.state_size, I.wstate);
+    return os;
+}
 // ------------------------------------------------------------
+// Class RK2Integrator
+// ------------------------------------------------------------
+// Default Constructor
+SA::RK2Integrator::RK2Integrator()
+    : FirstOrderODEVariableStepIntegrator() {
+    wstate = new double[state_size];
+    for (unsigned int i = 0; i < 2 ; i++) {
+        derivs[i] = new double[state_size];
+    }
+}
+// Constructor
 SA::RK2Integrator::RK2Integrator(
-    double h, int N, double* in_vars[], double* out_vars[], DerivsFunc dfunc, void* udata)
+    double h, unsigned int N, double* in_vars[], double* out_vars[], DerivsFunc dfunc, void* udata)
     : FirstOrderODEVariableStepIntegrator(h, N, in_vars, out_vars, dfunc, udata) {
     wstate = new double[N];
     for (unsigned int i = 0; i < 2 ; i++) {
         derivs[i] = new double[N];
     }
 }
+// Copy Constructor
+SA::RK2Integrator::RK2Integrator(const RK2Integrator& other)
+    : FirstOrderODEVariableStepIntegrator(other) {
+    wstate = new double[state_size];
+    std::copy(other.wstate, other.wstate+state_size, wstate);
+    for (unsigned int i=0 ; i<2; i++ ) {
+        derivs[i] = new double[state_size];
+        std::copy(other.derivs[i], other.derivs[i]+state_size, derivs[i]);
+    }
+}
+// Assignment Operator
+SA::RK2Integrator& SA::RK2Integrator::operator=( const SA::RK2Integrator& rhs) {
+    if (this != &rhs) {
+        // Call base assignment operator
+        FirstOrderODEVariableStepIntegrator::operator=(rhs);
+        // Duplicate rhs arrays
+        double* new_wstate = new double[rhs.state_size];
+        std::copy(rhs.wstate, rhs.wstate + rhs.state_size, new_wstate);
+        double* new_derivs[2];
+        for (unsigned int i=0 ; i<2; i++ ) {
+            new_derivs[i] = new double[rhs.state_size];
+            std::copy(rhs.derivs[i], rhs.derivs[i] + rhs.state_size, new_derivs[i]);
+        }
+        // Delete lhs arrays & replace with rhs arrays
+        if (wstate != NULL) delete wstate;
+        wstate = new_wstate;
+        for (unsigned int i=0 ; i<2; i++ ) {
+            if (derivs[i] != NULL) delete derivs[i];
+            derivs[i] = new_derivs[i];
+        }
+    }
+    return *this;
+}
+// Destructor
 SA::RK2Integrator::~RK2Integrator() {
-    delete (wstate);
+    if (wstate != NULL) delete (wstate);
     for (unsigned int i = 0; i < 2 ; i++) {
-        delete derivs[i];
+        if (derivs[i] != NULL) delete derivs[i];
     }
 }
 void SA::RK2Integrator::variable_step( double h) {
@@ -220,9 +531,36 @@ void SA::RK2Integrator::variable_step( double h) {
     }
     advanceIndyVar(h);
 }
+// Insertion Operator
+std::ostream& SA::operator<<(std::ostream& os, const RK2Integrator& I) {
+    os << (SA::FirstOrderODEVariableStepIntegrator)I;
+    os << "\n--- RK2Integrator ---";
+    os << "\nderivs    :[";
+    for (int i=0; i<2 ; i++) {
+        os << "\n";
+        stream_double_array(os, I.state_size, I.derivs[i]);
+    }
+    os << "\n]";
+    os << "\nwstate    :";
+    stream_double_array(os, I.state_size, I.wstate);
+    return os;
+}
 // ------------------------------------------------------------
+// Class RK4Integrator
+// ------------------------------------------------------------
+// Default Constructor
+SA::RK4Integrator::RK4Integrator()
+    : FirstOrderODEVariableStepIntegrator() {
+    for (unsigned int i = 0; i < 3 ; i++) {
+        wstate[i] = new double[state_size];
+    }
+    for (unsigned int i = 0; i < 4 ; i++) {
+        derivs[i] = new double[state_size];
+    }
+}
+// Constructor
 SA::RK4Integrator::RK4Integrator(
-    double h, int N, double* in_vars[], double* out_vars[], DerivsFunc dfunc, void* udata)
+    double h, unsigned int N, double* in_vars[], double* out_vars[], DerivsFunc dfunc, void* udata)
     : FirstOrderODEVariableStepIntegrator(h, N, in_vars, out_vars, dfunc, udata) {
     for (unsigned int i = 0; i < 3 ; i++) {
         wstate[i] = new double[N];
@@ -231,12 +569,53 @@ SA::RK4Integrator::RK4Integrator(
         derivs[i] = new double[N];
     }
 }
+// Copy Constructor
+SA::RK4Integrator::RK4Integrator(const RK4Integrator& other)
+    : FirstOrderODEVariableStepIntegrator(other) {
+    for (unsigned int i=0 ; i<3; i++ ) {
+        wstate[i] = new double[state_size];
+        std::copy(other.wstate[i], other.wstate[i]+state_size, wstate[i]);
+    }
+    for (unsigned int i=0 ; i<4; i++ ) {
+        derivs[i] = new double[state_size];
+        std::copy(other.derivs[i], other.derivs[i]+state_size, derivs[i]);
+    }
+}
+// Assignment Operator
+SA::RK4Integrator& SA::RK4Integrator::operator=( const SA::RK4Integrator& rhs) {
+    if (this != &rhs) {
+        // Call base assignment operator
+        FirstOrderODEVariableStepIntegrator::operator=(rhs);
+        // Duplicate rhs arrays
+        double* new_wstate[3];
+        for (unsigned int i=0 ; i<3; i++ ) {
+            new_wstate[i] = new double[rhs.state_size];
+            std::copy(rhs.wstate[i], rhs.wstate[i] + rhs.state_size, new_wstate[i]);
+        }
+        double* new_derivs[4];
+        for (unsigned int i=0 ; i<4; i++ ) {
+            new_derivs[i] = new double[rhs.state_size];
+            std::copy(rhs.derivs[i], rhs.derivs[i]+ rhs.state_size, new_derivs[i]);
+        }
+        // Delete lhs arrays & replace with rhs arrays
+        for (unsigned int i=0 ; i<3; i++ ) {
+            if (wstate[i] != NULL) delete wstate[i];
+            wstate[i] = new_wstate[i];
+        }
+        for (unsigned int i=0 ; i<4; i++ ) {
+            if (derivs[i] != NULL) delete derivs[i];
+            derivs[i] = new_derivs[i];
+        }
+    }
+    return *this;
+}
+// Destructor
 SA::RK4Integrator::~RK4Integrator() {
     for (unsigned int i = 0; i < 3 ; i++) {
-        delete (wstate[i]);
+        if (wstate[i] != NULL) delete (wstate[i]);
     }
     for (unsigned int i = 0; i < 4 ; i++) {
-        delete (derivs[i]);
+        if (derivs[i] != NULL) delete (derivs[i]);
     }
 }
 void SA::RK4Integrator::variable_step( double h) {
@@ -261,9 +640,40 @@ void SA::RK4Integrator::variable_step( double h) {
     }
     advanceIndyVar(h);
 }
+// Insertion Operator
+std::ostream& SA::operator<<(std::ostream& os, const RK4Integrator& I) {
+    os << (SA::FirstOrderODEVariableStepIntegrator)I;
+    os << "\n--- RK4Integrator ---";
+    os << "\nderivs    :[";
+    for (int i=0; i<4 ; i++) {
+        os << "\n";
+        stream_double_array(os, I.state_size, I.derivs[i]);
+    }
+    os << "\n]";
+    os << "\nwstate    :[";
+    for (int i=0; i<3 ; i++) {
+        os << "\n";
+        stream_double_array(os, I.state_size, I.wstate[i]);
+    }
+    os << "\n]";
+    return os;
+}
 // ------------------------------------------------------------
+// Class RK3_8Integrator
+// ------------------------------------------------------------
+// Default Constructor
+SA::RK3_8Integrator::RK3_8Integrator()
+    : FirstOrderODEVariableStepIntegrator() {
+    for (unsigned int i = 0; i < 3 ; i++) {
+        wstate[i] = new double[state_size];
+    }
+    for (unsigned int i = 0; i < 4 ; i++) {
+        derivs[i] = new double[state_size];
+    }
+}
+// Constructor
 SA::RK3_8Integrator::RK3_8Integrator(
-    double h, int N, double* in_vars[], double* out_vars[], DerivsFunc dfunc, void* udata)
+    double h, unsigned int N, double* in_vars[], double* out_vars[], DerivsFunc dfunc, void* udata)
     : FirstOrderODEVariableStepIntegrator(h, N, in_vars, out_vars ,dfunc, udata) {
     for (unsigned int i = 0; i < 3 ; i++) {
         wstate[i] = new double[N];
@@ -272,12 +682,53 @@ SA::RK3_8Integrator::RK3_8Integrator(
         derivs[i] = new double[N];
     }
 }
+// Copy Constructor
+SA::RK3_8Integrator::RK3_8Integrator(const RK3_8Integrator& other)
+    : FirstOrderODEVariableStepIntegrator(other) {
+    for (unsigned int i=0 ; i<3; i++ ) {
+        wstate[i] = new double[state_size];
+        std::copy(other.wstate[i], other.wstate[i]+state_size, wstate[i]);
+    }
+    for (unsigned int i=0 ; i<4; i++ ) {
+        derivs[i] = new double[state_size];
+        std::copy(other.derivs[i], other.derivs[i]+state_size, derivs[i]);
+    }
+}
+// Assignment Operator
+SA::RK3_8Integrator& SA::RK3_8Integrator::operator=( const SA::RK3_8Integrator& rhs) {
+    if (this != &rhs) {
+        // Call base assignment operator
+        FirstOrderODEVariableStepIntegrator::operator=(rhs);
+        // Duplicate rhs arrays
+        double* new_wstate[3];
+        for (unsigned int i=0 ; i<3; i++ ) {
+            new_wstate[i] = new double[rhs.state_size];
+            std::copy(rhs.wstate[i], rhs.wstate[i] + rhs.state_size, new_wstate[i]);
+        }
+        double* new_derivs[4];
+        for (unsigned int i=0 ; i<4; i++ ) {
+            new_derivs[i] = new double[rhs.state_size];
+            std::copy(rhs.derivs[i], rhs.derivs[i] + rhs.state_size, new_derivs[i]);
+        }
+        // Delete lhs arrays & replace with rhs arrays
+        for (unsigned int i=0 ; i<3; i++ ) {
+            if (wstate[i] != NULL) delete wstate[i];
+            wstate[i] = new_wstate[i];
+        }
+        for (unsigned int i=0 ; i<4; i++ ) {
+            if (derivs[i] != NULL) delete derivs[i];
+            derivs[i] = new_derivs[i];
+        }
+    }
+    return *this;
+}
+// Destructor
 SA::RK3_8Integrator::~RK3_8Integrator() {
     for (unsigned int i = 0; i < 3 ; i++) {
-        delete (wstate[i]);
+        if (wstate[i] != NULL) delete (wstate[i]);
     }
     for (unsigned int i = 0; i < 4 ; i++) {
-        delete (derivs[i]);
+        if (derivs[i] != NULL) delete (derivs[i]);
     }
 }
 void SA::RK3_8Integrator::variable_step( double h) {
@@ -302,73 +753,194 @@ void SA::RK3_8Integrator::variable_step( double h) {
     }
     advanceIndyVar(h);
 }
+// Insertion Operator
+std::ostream& SA::operator<<(std::ostream& os, const RK3_8Integrator& I) {
+    os << (SA::FirstOrderODEVariableStepIntegrator)I;
+    os << "\n--- RK3_8Integrator ---";
+    os << "\nderivs    :[";
+    for (int i=0; i<4 ; i++) {
+        os << "\n";
+        stream_double_array(os, I.state_size, I.derivs[i]);
+    }
+    os << "\n]";
+    os << "\nwstate    :[";
+    for (int i=0; i<3 ; i++) {
+        os << "\n";
+        stream_double_array(os, I.state_size, I.wstate[i]);
+    }
+    os << "\n]";
+    return os;
+}
+
 // ------------------------------------------------------------
+// Class EulerCromerIntegrator
+// ------------------------------------------------------------
+// Default Constructor
+SA::EulerCromerIntegrator::EulerCromerIntegrator()
+    :Integrator() {
+    nDimensions  = 1;
+    last_h = 0.0;
+    pos_p = NULL;
+    pos_in  = new double[nDimensions];
+    pos_out = new double[nDimensions];
+    vel_p = NULL;
+    vel_in  = new double[nDimensions];
+    vel_out = new double[nDimensions];
+    g_out   = new double[nDimensions];
+    gderivs = nil_gfunc;
+}
+// Constructor
 SA::EulerCromerIntegrator::EulerCromerIntegrator(
-    double dt, int N, double* xp[], double* vp[], DerivsFunc gfunc, DerivsFunc ffunc, void* udata)
+    double dt, unsigned int N, double* xp[], double* vp[], Derivs2Func gfunc, void* udata)
     :Integrator(dt, udata) {
     nDimensions  = N;
-    pos_p   = new double*[N];
+    last_h = 0.0;
+    pos_p = xp;
     pos_in  = new double[N];
     pos_out = new double[N];
-    vel_p   = new double*[N];
+    vel_p = vp;
     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++ ) {
-        pos_p[i] = xp[i];
-        vel_p[i] = vp[i];
-    }
     gderivs = gfunc;
-    fderivs = ffunc;
 }
+// Copy Constructor
+SA::EulerCromerIntegrator::EulerCromerIntegrator(const EulerCromerIntegrator& other)
+: Integrator(other.default_h, other.user_data) {
+    nDimensions = other.nDimensions;
+    last_h = other.last_h;
+    pos_p = other.pos_p;
+    pos_in = new double[nDimensions];
+    std::copy( other.pos_in, other.pos_in + nDimensions, pos_in);
+    pos_out = new double[nDimensions];
+    std::copy( other.pos_out, other.pos_out + nDimensions, pos_out);
+    vel_p = other.vel_p;
+    std::copy( other.vel_p, other.vel_p + nDimensions, vel_p);
+    vel_in = new double[nDimensions];
+    std::copy( other.vel_in, other.vel_in + nDimensions, vel_in);
+    vel_out = new double[nDimensions];
+    std::copy( other.vel_out, other.vel_out + nDimensions, vel_out);
+    g_out = new double[nDimensions];
+    std::copy( other.g_out, other.g_out + nDimensions, g_out);
+    gderivs = other.gderivs;
+}
+// Assignment Operator
+SA::EulerCromerIntegrator& SA::EulerCromerIntegrator::operator=( const SA::EulerCromerIntegrator& rhs) {
+    if (this != &rhs) {
+        // Call base assignment operator
+        Integrator::operator=(rhs);
+        // Duplicate the rhs arrays
+        double* new_pos_in  = new double[rhs.nDimensions];
+        std::copy( rhs.pos_in, rhs.pos_in + rhs.nDimensions, new_pos_in);
+        double* new_pos_out = new double[rhs.nDimensions];
+        std::copy( rhs.pos_out, rhs.pos_out + rhs.nDimensions, new_pos_out);
+        double* new_vel_in  = new double[rhs.nDimensions];
+        std::copy( rhs.vel_in, rhs.vel_in + rhs.nDimensions, new_vel_in);
+        double* new_vel_out = new double[rhs.nDimensions];
+        std::copy( rhs.vel_out, rhs.vel_out + rhs.nDimensions, new_vel_out);
+        double* new_g_out = new double[rhs.nDimensions];
+        std::copy( rhs.g_out, rhs.g_out + rhs.nDimensions, new_g_out);
+        // Delete lhs arrays & replace with rhs arrays
+        if (pos_in != NULL) delete pos_in;
+        pos_in = new_pos_in;
+        if (pos_out != NULL) delete pos_out;
+        pos_out = new_pos_out;
+        if (vel_in != NULL) delete vel_in;
+        vel_in = new_vel_in;
+        if (vel_out != NULL) delete vel_out;
+        vel_out = new_vel_out;
+        if (g_out != NULL) delete g_out;
+        g_out = new_g_out;
+        // Copy primitive members
+        last_h = rhs.last_h;
+        pos_p = rhs.pos_p;
+        vel_p = rhs.vel_p;
+        nDimensions = rhs.nDimensions;
+        gderivs = rhs.gderivs;
+    }
+    return *this;
+}
+// Destructor
 SA::EulerCromerIntegrator::~EulerCromerIntegrator() {
-
-    delete pos_in;
-    delete pos_out;
-    delete vel_in;
-    delete vel_out;
-    delete g_out;
-    delete f_out;
+    if (pos_in != NULL)  delete pos_in;
+    if (pos_out != NULL) delete pos_out;
+    if (vel_in != NULL)  delete vel_in;
+    if (vel_out != NULL) delete vel_out;
+    if (g_out != NULL)   delete g_out;
+}
+void SA::EulerCromerIntegrator::advanceIndyVar(double h) {
+    last_h = h; X_out = X_in + h;
 }
 void SA::EulerCromerIntegrator::step( double dt) {
-    (*gderivs)( X_in, pos_in, g_out, user_data);
+    (*gderivs)( X_in, pos_in, vel_in, g_out, user_data);
     for (int i=0 ; i<nDimensions; i++ ) {
         vel_out[i] = vel_in[i] + g_out[i] * dt;
+        pos_out[i] = pos_in[i] + vel_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;
+    advanceIndyVar(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<nDimensions; i++ ) {
-        pos_in[i] = *(pos_p[i]);
-        vel_in[i] = *(vel_p[i]);
+    if ((pos_in != NULL) && (vel_in != NULL)) {
+        for (int i=0 ; i<nDimensions; i++ ) {
+            pos_in[i] = *(pos_p[i]);
+            vel_in[i] = *(vel_p[i]);
+        }
+        Integrator::load();
+    } else {
+        std::cerr << "Error: SA::EulerCromerIntegrator::load(). Input variable pointers (pos_in, vel_in) aren't set." << std::endl;
     }
 }
 void SA::EulerCromerIntegrator::unload(){
-    for (int i=0 ; i<nDimensions; i++ ) {
-        *(pos_p[i]) = pos_out[i];
-        *(vel_p[i]) = vel_out[i];
+    if ((pos_in != NULL) && (vel_in != NULL)) {
+        for (int i=0 ; i<nDimensions; i++ ) {
+            *(pos_p[i]) = pos_out[i];
+            *(vel_p[i]) = vel_out[i];
+        }
+    } else {
+        std::cerr << "Error: SA::EulerCromerIntegrator::load(). Input variable pointers (pos_in, vel_in) aren't set." << std::endl;
     }
 }
 double SA::EulerCromerIntegrator::undo_integrate() {
-    for (int i=0 ; i<nDimensions; i++ ) {
-        *(pos_p[i]) = pos_in[i];
-        *(vel_p[i]) = vel_in[i];
+    if ((pos_in != NULL) && (vel_in != NULL)) {
+        for (unsigned int i=0 ; i<nDimensions; i++ ) {
+            *(pos_p[i]) = pos_in[i];
+            *(vel_p[i]) = vel_in[i];
+        }
+        std::copy(pos_in, pos_in + nDimensions, pos_out);
+        std::copy(vel_in, vel_in + nDimensions, vel_out);
+        return (SA::Integrator::undo_integrate());
+    } else {
+        std::cerr << "Error: SA::EulerCromerIntegrator::undo_integrate(). Input variable pointers (pos_in, vel_in) aren't set." << std::endl;
     }
-    return (SA::Integrator::undo_integrate());
+    return (last_h);
+}
+// Insertion Operator
+std::ostream& SA::operator<<(std::ostream& os, const EulerCromerIntegrator& I) {
+    os << (SA::Integrator)I;
+    os << "\n--- EulerCromerIntegrator ---";
+    os << "\nnDimensions : " << I.nDimensions;
+    os << "\nlast_h      : " << I.last_h;
+    os << "\npos_p       :" << I.pos_p;
+    os << "\nvel_p       :" << I.vel_p;
+    os << "\npos_in    :";
+    stream_double_array(os, I.nDimensions, I.pos_in);
+    os << "\npos_out    :";
+    stream_double_array(os, I.nDimensions, I.pos_out);
+    os << "\nvel_in    :";
+    stream_double_array(os, I.nDimensions, I.vel_in);
+    os << "\nvel_out    :";
+    stream_double_array(os, I.nDimensions, I.vel_out);
+    os << "\ng_out    :";
+    stream_double_array(os, I.nDimensions, I.g_out);
+    return os;
 }
 
+// ------------------------------------------------------------
+// Class ABMIntegrator
+// ------------------------------------------------------------
 // // Adams-Bashforth Coefficients
 static const double ABCoeffs[5][5] = {
   {1.0,           0.0,           0.0,          0.0,           0.0},
@@ -386,9 +958,15 @@ static const double AMCoeffs[5][5] = {
   {(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) {
 
+// Default Constructor
+
+// Constructor
+SA::ABMIntegrator::ABMIntegrator ( unsigned int history_size, double h, unsigned int N, double* in_vars[], double* out_vars[], DerivsFunc func, void* udata)
+    : FirstOrderODEIntegrator(h, N, in_vars, out_vars ,func, udata) {
+    if ((history_size < 1) || (history_size > 5)) {
+        throw std::invalid_argument("history_size must be in the range [1..5].");
+    }
     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;
@@ -399,6 +977,51 @@ SA::ABMIntegrator::ABMIntegrator ( int history_size, double h, int N, double* in
     }
     composite_deriv = new double[state_size];
 }
+// Copy Constructor
+SA::ABMIntegrator::ABMIntegrator(const ABMIntegrator& other)
+: FirstOrderODEIntegrator(other) {
+    bufferSize = other.bufferSize;
+    algorithmHistorySize = other.algorithmHistorySize;
+    hix = other.hix;
+    currentHistorySize = other.currentHistorySize;
+    deriv_history = new double*[bufferSize];
+    for (unsigned int i=0; i<bufferSize ; i++) {
+        deriv_history[i] = new double[state_size];
+        std::copy(other.deriv_history[i], other.deriv_history[i]+state_size, deriv_history[i]);
+    }
+    composite_deriv = new double[state_size];
+    std::copy(other.composite_deriv, other.composite_deriv+state_size, composite_deriv);
+}
+// Assignment Operator
+SA::ABMIntegrator& SA::ABMIntegrator::operator=( const SA::ABMIntegrator& rhs) {
+    if (this != &rhs) {
+        // Call base assignment operator
+        FirstOrderODEIntegrator::operator=(rhs);
+        // Duplicate rhs arrays
+        double** new_deriv_history = new double*[rhs.bufferSize];
+        for (unsigned int i=0 ; i<rhs.bufferSize; i++ ) {
+            new_deriv_history[i] = new double[rhs.state_size];
+            std::copy(rhs.deriv_history[i], rhs.deriv_history[i] + rhs.state_size, new_deriv_history[i]);
+        }
+        double* new_composite_deriv = new double[rhs.state_size];
+        std::copy(rhs.composite_deriv, rhs.composite_deriv+rhs.state_size, new_composite_deriv);
+        // Delete lhs arrays & replace with rhs arrays
+        for (unsigned int i=0 ; i<bufferSize; i++ ) {
+            delete deriv_history[i];
+        }
+        delete deriv_history;
+        deriv_history = new_deriv_history;
+        delete composite_deriv;
+        composite_deriv = new_composite_deriv;
+        // Copy primitive members
+        bufferSize = rhs.bufferSize;
+        algorithmHistorySize = rhs.algorithmHistorySize;
+        hix = rhs.hix;
+        currentHistorySize = rhs.currentHistorySize;
+    }
+    return *this;
+}
+// Destructor
 SA::ABMIntegrator::~ABMIntegrator() {
     for (int i=0; i<bufferSize ; i++) {
         delete (deriv_history[i]);
@@ -433,20 +1056,39 @@ void SA::ABMIntegrator::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::ABMIntegrator::undo_integrate() {
     hix = (hix+bufferSize-1)%bufferSize;
     return (FirstOrderODEIntegrator::undo_integrate());
 }
+// Insertion Operator
+std::ostream& SA::operator<<(std::ostream& os, const ABMIntegrator& I) {
+    os << (SA::FirstOrderODEIntegrator)I;
+    os << "\n--- ABMIntegrator ---";
+    os << "\nbufferSize           : " << I.bufferSize;
+    os << "\nalgorithmHistorySize : " << I.algorithmHistorySize;
+    os << "\nhix                  : " << I.hix;
+    os << "\ncurrentHistorySize   : " << I.currentHistorySize;
+    os << "\nderiv_history        :[";
+    for (int i=0; i<I.bufferSize ; i++) {
+        os << "\n";
+        stream_double_array(os, I.state_size, I.deriv_history[i]);
+    }
+    os << "\n]";
+    os << "\ncomposite_deriv     :";
+    stream_double_array(os, I.state_size, I.composite_deriv);
+    return os;
+}
 
-// =======================================================================
+// ------------------------------------------------------------
+// Class ABM2Integrator
+// ------------------------------------------------------------
 static const double AB2Coeffs[2] = {(3/2.0), (-1/2.0)};
 static const double AM2Coeffs[2] = {(1/2.0),  (1/2.0)};
 
 SA::ABM2Integrator::ABM2Integrator (
-    double h, int N, double* in_vars[], double* out_vars[], DerivsFunc dfunc, void* udata)
+    double h, unsigned int N, double* in_vars[], double* out_vars[], DerivsFunc dfunc, void* udata)
     : FirstOrderODEIntegrator(h, N, in_vars, out_vars, dfunc, udata) {
     // The amount of history that we intend to use in this ABM integrator.
     bufferSize=2+1; // The size of the buffer needs to be one more than the history so that an integration step can be reset.
@@ -469,6 +1111,52 @@ SA::ABM2Integrator::ABM2Integrator (
     }
     priming_integrator = new SA::RK2Integrator(h, N, priming_integrator_in_vars, priming_integrator_out_vars, dfunc, udata);
 }
+
+
+SA::ABM2Integrator::ABM2Integrator(const ABM2Integrator& other)
+: FirstOrderODEIntegrator(other) {
+    bufferSize = other.bufferSize;
+    hix = other.hix;
+    currentHistorySize = other.currentHistorySize;
+    deriv_history = new double*[bufferSize];
+    for (unsigned int i=0; i<bufferSize ; i++) {
+        deriv_history[i] = new double[state_size];
+        std::copy(other.deriv_history[i], other.deriv_history[i]+state_size, deriv_history[i]);
+    }
+    composite_deriv = new double[state_size];
+    std::copy(other.composite_deriv, other.composite_deriv+state_size, composite_deriv);
+
+    priming_integrator = new RK2Integrator(*other.priming_integrator);
+}
+// Assignment Operator
+SA::ABM2Integrator& SA::ABM2Integrator::operator=( const SA::ABM2Integrator& rhs) {
+    if (this != &rhs) {
+        // Call base assignment operator
+        FirstOrderODEIntegrator::operator=(rhs);
+        // Duplicate rhs arrays
+        double** new_deriv_history = new double*[rhs.bufferSize];
+        for (unsigned int i=0 ; i<rhs.bufferSize; i++ ) {
+            new_deriv_history[i] = new double[rhs.state_size];
+            std::copy(rhs.deriv_history[i], rhs.deriv_history[i] + rhs.state_size, new_deriv_history[i]);
+        }
+        double* new_composite_deriv = new double[rhs.state_size];
+        std::copy(rhs.composite_deriv, rhs.composite_deriv+rhs.state_size, new_composite_deriv);
+        // Delete lhs arrays & replace with rhs arrays
+        for (unsigned int i=0 ; i<bufferSize; i++ ) {
+            delete deriv_history[i];
+        }
+        delete deriv_history;
+        deriv_history = new_deriv_history;
+        delete composite_deriv;
+        composite_deriv = new_composite_deriv;
+        // Copy primitive members
+        bufferSize = rhs.bufferSize;
+        hix = rhs.hix;
+        currentHistorySize = rhs.currentHistorySize;
+    }
+    return *this;
+}
+
 SA::ABM2Integrator::~ABM2Integrator() {
     for (int i=0; i<bufferSize ; i++) {
         delete (deriv_history[i]);
@@ -518,12 +1206,32 @@ double SA::ABM2Integrator::undo_integrate() {
     hix = (hix+bufferSize-1)%bufferSize;
     return ( FirstOrderODEIntegrator::undo_integrate());
 }
-// =======================================================================
+
+// Insertion Operator
+std::ostream& SA::operator<<(std::ostream& os, const ABM2Integrator& I) {
+    os << (SA::FirstOrderODEIntegrator)I;
+    os << "\n--- ABM2Integrator ---";
+    os << "\nbufferSize           : " << I.bufferSize;
+    os << "\nhix                  : " << I.hix;
+    os << "\ncurrentHistorySize   : " << I.currentHistorySize;
+    os << "\nderiv_history        :[";
+    for (int i=0; i<I.bufferSize ; i++) {
+        os << "\n";
+        stream_double_array(os, I.state_size, I.deriv_history[i]);
+    }
+    os << "\n]";
+    os << "\ncomposite_deriv     :";
+    stream_double_array(os, I.state_size, I.composite_deriv);
+    return os;
+}
+// ------------------------------------------------------------
+// Class ABM4Integrator
+// ------------------------------------------------------------
 static const double AB4Coeffs[4] = {(55/24.0), (-59/24.0), (37/24.0), (-9/24.0)};
 static const double AM4Coeffs[4] = {(9/24.0), (19/24.0), (-5/24.0), (1/24.0)};
 
 SA::ABM4Integrator::ABM4Integrator (
-    double h, int N, double* in_vars[], double* out_vars[], DerivsFunc dfunc, void* udata)
+    double h, unsigned int N, double* in_vars[], double* out_vars[], DerivsFunc dfunc, void* udata)
     : FirstOrderODEIntegrator(h, N, in_vars, out_vars, dfunc, udata) {
 
     // The amount of history that we intend to use in this ABM integrator.
@@ -546,6 +1254,50 @@ SA::ABM4Integrator::ABM4Integrator (
     }
     priming_integrator = new SA::RK4Integrator(h, N, priming_integrator_in_vars, priming_integrator_out_vars, dfunc, udata);
 }
+
+SA::ABM4Integrator::ABM4Integrator(const ABM4Integrator& other)
+: FirstOrderODEIntegrator(other) {
+    bufferSize = other.bufferSize;
+    hix = other.hix;
+    currentHistorySize = other.currentHistorySize;
+    deriv_history = new double*[bufferSize];
+    for (unsigned int i=0; i<bufferSize ; i++) {
+        deriv_history[i] = new double[state_size];
+        std::copy(other.deriv_history[i], other.deriv_history[i]+state_size, deriv_history[i]);
+    }
+    composite_deriv = new double[state_size];
+    std::copy(other.composite_deriv, other.composite_deriv+state_size, composite_deriv);
+
+    priming_integrator = new RK4Integrator(*other.priming_integrator);
+}
+// Assignment Operator
+SA::ABM4Integrator& SA::ABM4Integrator::operator=( const SA::ABM4Integrator& rhs) {
+    if (this != &rhs) {
+        // Call base assignment operator
+        FirstOrderODEIntegrator::operator=(rhs);
+        // Duplicate rhs arrays
+        double** new_deriv_history = new double*[rhs.bufferSize];
+        for (unsigned int i=0 ; i<rhs.bufferSize; i++ ) {
+            new_deriv_history[i] = new double[rhs.state_size];
+            std::copy(rhs.deriv_history[i], rhs.deriv_history[i] + rhs.state_size, new_deriv_history[i]);
+        }
+        double* new_composite_deriv = new double[rhs.state_size];
+        std::copy(rhs.composite_deriv, rhs.composite_deriv+rhs.state_size, new_composite_deriv);
+        // Delete lhs arrays & replace with rhs arrays
+        for (unsigned int i=0 ; i<bufferSize; i++ ) {
+            delete deriv_history[i];
+        }
+        delete deriv_history;
+        deriv_history = new_deriv_history;
+        delete composite_deriv;
+        composite_deriv = new_composite_deriv;
+        // Copy primitive members
+        bufferSize = rhs.bufferSize;
+        hix = rhs.hix;
+        currentHistorySize = rhs.currentHistorySize;
+    }
+    return *this;
+}
 SA::ABM4Integrator::~ABM4Integrator() {
     for (int i=0; i<bufferSize ; i++) {
         delete (deriv_history[i]);
@@ -594,3 +1346,20 @@ double SA::ABM4Integrator::undo_integrate() {
     hix = (hix+bufferSize-1)%bufferSize;
     return (FirstOrderODEIntegrator::undo_integrate());
 }
+// Insertion Operator
+std::ostream& SA::operator<<(std::ostream& os, const ABM4Integrator& I) {
+    os << (SA::FirstOrderODEIntegrator)I;
+    os << "\n--- ABM4Integrator ---";
+    os << "\nbufferSize           : " << I.bufferSize;
+    os << "\nhix                  : " << I.hix;
+    os << "\ncurrentHistorySize   : " << I.currentHistorySize;
+    os << "\nderiv_history        :[";
+    for (int i=0; i<I.bufferSize ; i++) {
+        os << "\n";
+        stream_double_array(os, I.state_size, I.deriv_history[i]);
+    }
+    os << "\n]";
+    os << "\ncomposite_deriv     :";
+    stream_double_array(os, I.state_size, I.composite_deriv);
+    return os;
+}
diff --git a/trick_source/trick_utils/SAIntegrator/unittest/EulerIntegrator_unittest.cc b/trick_source/trick_utils/SAIntegrator/unittest/EulerIntegrator_unittest.cc
new file mode 100644
index 00000000..cb96ea9f
--- /dev/null
+++ b/trick_source/trick_utils/SAIntegrator/unittest/EulerIntegrator_unittest.cc
@@ -0,0 +1,86 @@
+#include <gtest/gtest.h>
+#include <iostream>
+#include "SAIntegrator.hh"
+#include <math.h>
+
+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(EulerIntegrator_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, 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/FirstOrderODEIntegrator_unittest.cc b/trick_source/trick_utils/SAIntegrator/unittest/FirstOrderODEIntegrator_unittest.cc
index 022358f8..e617b80d 100644
--- a/trick_source/trick_utils/SAIntegrator/unittest/FirstOrderODEIntegrator_unittest.cc
+++ b/trick_source/trick_utils/SAIntegrator/unittest/FirstOrderODEIntegrator_unittest.cc
@@ -156,4 +156,65 @@ TEST(FirstOrderODEIntegrator_unittest, test_4) {
 
     double x = integ.getIndyVar();
     EXPECT_NEAR(x, 0.6, EXCEPTABLE_ERROR);
+
+}
+
+TEST(FirstOrderODEIntegrator_unittest, test_copy_constructor) {
+    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 integ1( h, N, varptrs, varptrs, calc_derivs, NULL);
+
+    // Do an integration step to change array values to non-zero.
+    integ1.integrate();
+
+    // Use the copy-constructor to create create a duplicate FirstOrderODEIntegrator.
+    SA::FirstOrderODEIntegrator integ2(integ1);
+
+    // We don't want to unload integ1 here because it will then effect integ2,
+    // since both are pointing to the same input/output arrays.
+    integ1.load(); integ1.step();
+    integ2.load(); integ2.step();
+
+    // Create a text representation for each of the FirstOrderODEIntegrators.
+    std::stringstream ss1;
+    std::stringstream ss2;
+    ss1 << integ1 << std::endl;
+    ss2 << integ2 << std::endl;
+
+    // Compare the text representations. They should be identical.
+    int result = ss1.str().compare(ss2.str());
+    EXPECT_EQ(result, 0);
+}
+
+TEST(FirstOrderODEIntegrator_unittest, test_assignment_operator) {
+    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 integ1( h, N, varptrs, varptrs, calc_derivs, NULL);
+
+    // Do an integration step to change array values to non-zero.
+    integ1.integrate();
+
+    // Use the assignment-operator to create a duplicate FirstOrderODEIntegrator.
+    SA::FirstOrderODEIntegrator integ2 = integ1;
+
+    // We don't want to unload integ1 here because it will then effect integ2,
+    // since both are pointing to the same input/output arrays.
+    integ1.load(); integ1.step();
+    integ2.load(); integ2.step();
+
+    // Create a text representation for each of the FirstOrderODEIntegrators.
+    std::stringstream ss1;
+    std::stringstream ss2;
+    ss1 << integ1 << std::endl;
+    ss2 << integ2 << std::endl;
+
+    // Compare the text representations. They should be identical.
+    int result = ss1.str().compare(ss2.str());
+    EXPECT_EQ(result, 0);
 }
diff --git a/trick_source/trick_utils/SAIntegrator/unittest/FirstOrderODEVariableStepIntegrator_unittest.cc b/trick_source/trick_utils/SAIntegrator/unittest/FirstOrderODEVariableStepIntegrator_unittest.cc
index 14675eac..5f5b56e1 100644
--- a/trick_source/trick_utils/SAIntegrator/unittest/FirstOrderODEVariableStepIntegrator_unittest.cc
+++ b/trick_source/trick_utils/SAIntegrator/unittest/FirstOrderODEVariableStepIntegrator_unittest.cc
@@ -38,69 +38,6 @@ void deriv1( double t __attribute__((unused)),
 
 #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};
@@ -108,11 +45,11 @@ TEST(FirstOrderODEVariableStepIntegrator_unittest, RungeKutta38_1) {
    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);
+   SA::RK3_8Integrator integ1(dt, 4, state_var_p, state_var_p, deriv4, NULL);
    while (t < 2.0) {
-        integ.load();
-        integ.step();
-        integ.unload();
+        integ1.load();
+        integ1.step();
+        integ1.unload();
         count ++;
         t = count * dt;
     }
@@ -120,53 +57,17 @@ TEST(FirstOrderODEVariableStepIntegrator_unittest, RungeKutta38_1) {
     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);
+
+    // Use the copy-constructor to create create a duplicate RK3_8Integrator.
+    SA::RK3_8Integrator integ2(integ1);
+
+    // Create a text representation for each of the RK3_8Integrator.
+    std::stringstream ss1;
+    std::stringstream ss2;
+    ss1 << integ1 << std::endl;
+    ss2 << integ2 << std::endl;
+
+    // Compare the text representations. They should be identical.
+    int result = ss1.str().compare(ss2.str());
+    EXPECT_EQ(result, 0);
 }
diff --git a/trick_source/trick_utils/SAIntegrator/unittest/Makefile b/trick_source/trick_utils/SAIntegrator/unittest/Makefile
index 1a078d94..c67a288a 100644
--- a/trick_source/trick_utils/SAIntegrator/unittest/Makefile
+++ b/trick_source/trick_utils/SAIntegrator/unittest/Makefile
@@ -13,12 +13,25 @@ SAI_LIBOBJS = ${SAI_OBJDIR}/Integrator.o
 
 LIBDIRS = -L${SAI_LIBDIR} -L${GTEST_HOME}/lib64 -L${GTEST_HOME}/lib
 
+TEST_PROGRAMS = SAIntegrator_unittest\
+	            FirstOrderODEIntegrator_unittest\
+	            FirstOrderODEVariableStepIntegrator_unittest\
+	            EulerIntegrator_unittest\
+				RK2Integrator_unittest\
+				RK4Integrator_unittest\
+				RK3_8Integrator_unittest\
+	            RootFinder_unittest
+
 all: test
 
-test: SAIntegrator_unittest FirstOrderODEIntegrator_unittest FirstOrderODEVariableStepIntegrator_unittest RootFinder_unittest
+test: $(TEST_PROGRAMS)
 	./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
+	./EulerIntegrator_unittest --gtest_output=xml:${TRICK_HOME}/trick_test/EulerIntegrator_unittest.xml
+	./RK2Integrator_unittest --gtest_output=xml:${TRICK_HOME}/trick_test/RK2Integrator_unittest.xml
+	./RK4Integrator_unittest --gtest_output=xml:${TRICK_HOME}/trick_test/RK4Integrator_unittest.xml
+	./RK3_8Integrator_unittest --gtest_output=xml:${TRICK_HOME}/trick_test/RK3_8Integrator_unittest.xml
 	./RootFinder_unittest --gtest_output=xml:${TRICK_HOME}/trick_test/RootFinder_unittest.xml
 
 SAIntegrator_unittest.o : SAIntegrator_unittest.cc
@@ -38,7 +51,31 @@ FirstOrderODEVariableStepIntegrator_unittest.o : FirstOrderODEVariableStepIntegr
 
 FirstOrderODEVariableStepIntegrator_unittest : ${SAI_LIBDIR}/${SAI_LIBNAME} FirstOrderODEVariableStepIntegrator_unittest.o
 	$(TRICK_CXX) $(TRICK_CPPFLAGS) -o $@ $^ ${LIBDIRS} -lSAInteg -lgtest -lgtest_main -lpthread
+# ====
+EulerIntegrator_unittest.o : EulerIntegrator_unittest.cc
+	$(TRICK_CXX) $(TRICK_CPPFLAGS) $(INCLUDE_DIRS) -c $<
 
+EulerIntegrator_unittest : ${SAI_LIBDIR}/${SAI_LIBNAME} EulerIntegrator_unittest.o
+	$(TRICK_CXX) $(TRICK_CPPFLAGS) -o $@ $^ ${LIBDIRS} -lSAInteg -lgtest -lgtest_main -lpthread
+# ====
+RK2Integrator_unittest.o : RK2Integrator_unittest.cc
+	$(TRICK_CXX) $(TRICK_CPPFLAGS) $(INCLUDE_DIRS) -c $<
+
+RK2Integrator_unittest : ${SAI_LIBDIR}/${SAI_LIBNAME} RK2Integrator_unittest.o
+	$(TRICK_CXX) $(TRICK_CPPFLAGS) -o $@ $^ ${LIBDIRS} -lSAInteg -lgtest -lgtest_main -lpthread
+# ====
+RK4Integrator_unittest.o : RK4Integrator_unittest.cc
+	$(TRICK_CXX) $(TRICK_CPPFLAGS) $(INCLUDE_DIRS) -c $<
+
+RK4Integrator_unittest : ${SAI_LIBDIR}/${SAI_LIBNAME} RK4Integrator_unittest.o
+	$(TRICK_CXX) $(TRICK_CPPFLAGS) -o $@ $^ ${LIBDIRS} -lSAInteg -lgtest -lgtest_main -lpthread
+# ====
+RK3_8Integrator_unittest.o : RK3_8Integrator_unittest.cc
+	$(TRICK_CXX) $(TRICK_CPPFLAGS) $(INCLUDE_DIRS) -c $<
+
+RK3_8Integrator_unittest : ${SAI_LIBDIR}/${SAI_LIBNAME} RK3_8Integrator_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 $<
 
@@ -52,4 +89,11 @@ clean:
 	${RM} *.o
 
 spotless: clean
+	${RM} SAIntegrator_unittest
+	${RM} FirstOrderODEIntegrator_unittest
 	${RM} FirstOrderODEVariableStepIntegrator_unittest
+	${RM} EulerIntegrator_unittest
+	${RM} RK2Integrator_unittest
+	${RM} RK4Integrator_unittest
+	${RM} RK3_8Integrator_unittest
+	${RM} RootFinder_unittest
diff --git a/trick_source/trick_utils/SAIntegrator/unittest/RK2Integrator_unittest.cc b/trick_source/trick_utils/SAIntegrator/unittest/RK2Integrator_unittest.cc
new file mode 100644
index 00000000..604589fd
--- /dev/null
+++ b/trick_source/trick_utils/SAIntegrator/unittest/RK2Integrator_unittest.cc
@@ -0,0 +1,116 @@
+#include <gtest/gtest.h>
+#include <iostream>
+#include "SAIntegrator.hh"
+#include <math.h>
+
+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;
+}
+
+#define EXCEPTABLE_ERROR 0.00000000001
+
+TEST(RK2Integrator_unittest, test_load_step_unload) {
+
+   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 integ1(dt, 4, state_var_p, state_var_p, deriv2, NULL);
+   while (t < 2.0) {
+        integ1.load();
+        integ1.step();
+        integ1.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(RK2Integrator_unittest, test_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 dt = 0.01;
+   unsigned int count = 0;
+   double t = count * dt;
+   SA::RK2Integrator integ1(dt, 4, state_var_p, state_var_p, deriv2, NULL);
+   while (t < 2.0) {
+        integ1.integrate();
+        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(RK2Integrator_unittest, copy_constructor) {
+
+   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 integ1(dt, 4, state_var_p, state_var_p, deriv2, NULL);
+   while (t < 2.0) {
+        integ1.integrate();
+        count ++;
+        t = count * dt;
+    }
+
+    // Use the copy-constructor to create create a duplicate RK2Integrator.
+    SA::RK2Integrator integ2(integ1);
+
+    // Create a text representation for each of the RK2Integrators.
+    std::stringstream ss1;
+    std::stringstream ss2;
+    ss1 << integ1 << std::endl;
+    ss2 << integ2 << std::endl;
+
+    // Compare the text representations. They should be identical.
+    int result = ss1.str().compare(ss2.str());
+    EXPECT_EQ(result, 0);
+}
+
+TEST(RK2Integrator_unittest, assignment_operator) {
+
+   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 integ1(dt, 4, state_var_p, state_var_p, deriv2, NULL);
+   while (t < 2.0) {
+        integ1.integrate();
+        count ++;
+        t = count * dt;
+    }
+
+    // Use the assignment operator to create create a duplicate RK2Integrator.
+    SA::RK2Integrator integ2;
+    integ2 = integ1;
+
+    // Create a text representation for each of the RK2Integrators.
+    std::stringstream ss1;
+    std::stringstream ss2;
+    ss1 << integ1 << std::endl;
+    ss2 << integ2 << std::endl;
+
+    // Compare the text representations. They should be identical.
+    int result = ss1.str().compare(ss2.str());
+    EXPECT_EQ(result, 0);
+}
diff --git a/trick_source/trick_utils/SAIntegrator/unittest/RK3_8Integrator_unittest.cc b/trick_source/trick_utils/SAIntegrator/unittest/RK3_8Integrator_unittest.cc
new file mode 100644
index 00000000..a2e6e433
--- /dev/null
+++ b/trick_source/trick_utils/SAIntegrator/unittest/RK3_8Integrator_unittest.cc
@@ -0,0 +1,52 @@
+#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;
+}
+
+#define EXCEPTABLE_ERROR 0.00000000001
+
+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 integ1(dt, 4, state_var_p, state_var_p, deriv4, NULL);
+   while (t < 2.0) {
+        integ1.load();
+        integ1.step();
+        integ1.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);
+
+    // Use the copy-constructor to create create a duplicate RK3_8Integrator.
+    SA::RK3_8Integrator integ2(integ1);
+
+    // Create a text representation for each of the RK3_8Integrator.
+    std::stringstream ss1;
+    std::stringstream ss2;
+    ss1 << integ1 << std::endl;
+    ss2 << integ2 << std::endl;
+
+    // Compare the text representations. They should be identical.
+    int result = ss1.str().compare(ss2.str());
+    EXPECT_EQ(result, 0);
+}
diff --git a/trick_source/trick_utils/SAIntegrator/unittest/RK4Integrator_unittest.cc b/trick_source/trick_utils/SAIntegrator/unittest/RK4Integrator_unittest.cc
new file mode 100644
index 00000000..9c319282
--- /dev/null
+++ b/trick_source/trick_utils/SAIntegrator/unittest/RK4Integrator_unittest.cc
@@ -0,0 +1,52 @@
+#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;
+}
+
+#define EXCEPTABLE_ERROR 0.00000000001
+
+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 integ1(dt, 4, state_var_p, state_var_p, deriv4, NULL);
+   while (t < 2.0) {
+        integ1.load();
+        integ1.step();
+        integ1.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);
+
+    // Use the copy-constructor to create create a duplicate RK4Integrator.
+    SA::RK4Integrator integ2(integ1);
+
+    // Create a text representation for each of the RK4Integrator.
+    std::stringstream ss1;
+    std::stringstream ss2;
+    ss1 << integ1 << std::endl;
+    ss2 << integ2 << std::endl;
+
+    // Compare the text representations. They should be identical.
+    int result = ss1.str().compare(ss2.str());
+    EXPECT_EQ(result, 0);
+}
diff --git a/trick_source/trick_utils/SAIntegrator/unittest/RootFinder_unittest.cc b/trick_source/trick_utils/SAIntegrator/unittest/RootFinder_unittest.cc
index 7ac35dbe..bb0cd33a 100644
--- a/trick_source/trick_utils/SAIntegrator/unittest/RootFinder_unittest.cc
+++ b/trick_source/trick_utils/SAIntegrator/unittest/RootFinder_unittest.cc
@@ -78,3 +78,12 @@ TEST(RootFinder_unittest, Negative_constraint_with_positive_slope_function) {
     root_error = root_finder->find_roots(9, func_positive_slope(9));
     EXPECT_NEAR(root_error, DBL_MAX, EXCEPTABLE_ERROR);
 }
+
+TEST(RootFinder_unittest, serialize) {
+    RootFinder* root_finder = new RootFinder ( 0.00000000001, Negative );
+    double root_error;
+    root_error = root_finder->find_roots(8, func_positive_slope(8));
+    EXPECT_NEAR(root_error, DBL_MAX, EXCEPTABLE_ERROR);
+    root_error = root_finder->find_roots(9, func_positive_slope(9));
+    EXPECT_NEAR(root_error, DBL_MAX, EXCEPTABLE_ERROR);
+}
diff --git a/trick_source/trick_utils/SAIntegrator/unittest/SAIntegrator_unittest.cc b/trick_source/trick_utils/SAIntegrator/unittest/SAIntegrator_unittest.cc
index 96c2a0d9..0088a56f 100644
--- a/trick_source/trick_utils/SAIntegrator/unittest/SAIntegrator_unittest.cc
+++ b/trick_source/trick_utils/SAIntegrator/unittest/SAIntegrator_unittest.cc
@@ -72,3 +72,35 @@ TEST(Integrator, test_undo_integrate) {
     x = integ.getIndyVar();
     EXPECT_NEAR(x, 4.0, EXCEPTABLE_ERROR);
 }
+
+TEST(Integrator, test_copy_constructor) {
+    std::stringstream ss1;
+    std::stringstream ss2;
+    double h = 0.1;
+    void* udata = NULL;
+    SA::Integrator  integ( h, udata);
+    integ.setIndyVar(4.0);
+    integ.integrate();
+    // Test Copy Constructor
+    SA::Integrator  integ2(integ);
+    ss1 << integ << std::endl;
+    ss2 << integ2 << std::endl;
+    int result = ss1.str().compare(ss2.str());
+    EXPECT_EQ(result, 0);
+}
+
+TEST(Integrator, test_assignment_operator) {
+    std::stringstream ss1;
+    std::stringstream ss2;
+    double h = 0.1;
+    void* udata = NULL;
+    SA::Integrator  integ( h, udata);
+    integ.setIndyVar(4.0);
+    integ.integrate();
+    // Test Copy Constructor
+    SA::Integrator integ2 = integ;
+    ss1 << integ << std::endl;
+    ss2 << integ2 << std::endl;
+    int result = ss1.str().compare(ss2.str());
+    EXPECT_EQ(result, 0);
+}