diff --git a/examples/Experimental/DynamicConstrainedOpt/DynamicConstrainedOpt.cpp b/examples/Experimental/DynamicConstrainedOpt/DynamicConstrainedOpt.cpp
index b8def3c1..9c3654b2 100644
--- a/examples/Experimental/DynamicConstrainedOpt/DynamicConstrainedOpt.cpp
+++ b/examples/Experimental/DynamicConstrainedOpt/DynamicConstrainedOpt.cpp
@@ -1,6 +1,6 @@
-
 #include <iomanip>
 #include <iostream>
+#include <memory>
 
 #include <IpIpoptApplication.hpp>
 #include <IpSolveStatistics.hpp>
@@ -20,47 +20,47 @@ int main()
   using namespace GridKit::Testing;
 
   // Create an infinite bus
-  BaseBus<double, size_t>* bus = new BusSlack<double, size_t>(1.0, 0.0);
+  BusSlack<double, size_t> bus(1.0, 0.0);
 
   // Attach a generator to that bus
-  Generator2<double, size_t>* gen = new Generator2<double, size_t>(bus);
+  Generator2<double, size_t> gen(&bus);
 
   // Create a system model
-  SystemModel<double, size_t>* model = new SystemModel<double, size_t>();
-  model->addBus(bus);
-  model->addComponent(gen);
+  SystemModel<double, size_t> model;
+  model.addBus(&bus);
+  model.addComponent(&gen);
 
   // allocate model components
-  model->allocate();
+  model.allocate();
 
   // Create numerical integrator and configure it for the generator model
-  Ida<double, size_t>* idas = new Ida<double, size_t>(model);
+  Ida<double, size_t> idas(&model);
 
   double t_init  = 0.0;
   double t_final = 20.0;
 
   // setup simulation
-  idas->configureSimulation();
-  idas->configureAdjoint();
-  idas->getDefaultInitialCondition();
-  idas->initializeSimulation(t_init);
-  idas->configureQuadrature();
-  idas->initializeQuadrature();
+  idas.configureSimulation();
+  idas.configureAdjoint();
+  idas.getDefaultInitialCondition();
+  idas.initializeSimulation(t_init);
+  idas.configureQuadrature();
+  idas.initializeQuadrature();
 
   double t_fault = 0.02;
   double t_clear = 0.06;
-  idas->runSimulation(t_fault);
+  idas.runSimulation(t_fault);
   // create initial condition after a fault
   {
-    gen->V() = 0.0;
-    idas->runSimulation(t_clear, 2);
-    gen->V()     = 1.0;
-    gen->theta() = -0.01;
-    idas->saveInitialCondition();
+    gen.V() = 0.0;
+    idas.runSimulation(t_clear, 2);
+    gen.V()     = 1.0;
+    gen.theta() = -0.01;
+    idas.saveInitialCondition();
   }
 
   // Set integration time for dynamic constrained optimization
-  idas->setIntegrationTime(t_init, t_final, 100);
+  idas.setIntegrationTime(t_init, t_final, 100);
 
   // Guess optimization parameter value
   double Pm = 0.7;
@@ -87,10 +87,10 @@ int main()
 
   // Create dynamic objective interface to Ipopt solver
   Ipopt::SmartPtr<Ipopt::TNLP> ipoptDynamicObjectiveInterface =
-      new IpoptInterface::DynamicObjective<double, size_t>(idas);
+      new IpoptInterface::DynamicObjective<double, size_t>(&idas);
 
   // Initialize problem
-  model->param()[0] = Pm;
+  model.param()[0] = Pm;
 
   // Solve the problem
   status = ipoptApp->OptimizeTNLP(ipoptDynamicObjectiveInterface);
@@ -101,24 +101,24 @@ int main()
     // Print result
     std::cout << "\nSucess:\n The problem solved in "
               << ipoptApp->Statistics()->IterationCount() << " iterations!\n"
-              << " Optimal value of Pm = " << model->param()[0] << "\n"
+              << " Optimal value of Pm = " << model.param()[0] << "\n"
               << " The final value of the objective function G(Pm) = "
               << ipoptApp->Statistics()->FinalObjective() << "\n\n";
   }
 
   // Store dynamic objective optimization results
-  double* results = new double[model->sizeParams()];
-  for (unsigned i = 0; i < model->sizeParams(); ++i)
+  double* results = new double[model.sizeParams()];
+  for (unsigned i = 0; i < model.sizeParams(); ++i)
   {
-    results[i] = model->param()[i];
+    results[i] = model.param()[i];
   }
 
   // Create dynamic constraint interface to Ipopt solver
   Ipopt::SmartPtr<Ipopt::TNLP> ipoptDynamicConstraintInterface =
-      new IpoptInterface::DynamicConstraint<double, size_t>(idas);
+      new IpoptInterface::DynamicConstraint<double, size_t>(&idas);
 
   // Initialize problem
-  model->param()[0] = Pm;
+  model.param()[0] = Pm;
 
   // Solve the problem
   status = ipoptApp->OptimizeTNLP(ipoptDynamicConstraintInterface);
@@ -129,16 +129,16 @@ int main()
     // Print result
     std::cout << "\nSucess:\n The problem solved in "
               << ipoptApp->Statistics()->IterationCount() << " iterations!\n"
-              << " Optimal value of Pm = " << model->param()[0] << "\n"
+              << " Optimal value of Pm = " << model.param()[0] << "\n"
               << " The final value of the objective function G(Pm) = "
               << ipoptApp->Statistics()->FinalObjective() << "\n\n";
   }
 
   // Compare results of the two optimization methods
   int retval = 0;
-  for (unsigned i = 0; i < model->sizeParams(); ++i)
+  for (unsigned i = 0; i < model.sizeParams(); ++i)
   {
-    if (!isEqual(results[i], model->param()[i], 10 * tol))
+    if (!isEqual(results[i], model.param()[i], 10 * tol))
       --retval;
   }
 
@@ -148,7 +148,5 @@ int main()
   }
 
   delete[] results;
-  delete idas;
-  delete model;
   return retval;
 }
diff --git a/examples/Experimental/GenInfiniteBus/GenInfiniteBus.cpp b/examples/Experimental/GenInfiniteBus/GenInfiniteBus.cpp
index b3ea77a0..143e3d09 100644
--- a/examples/Experimental/GenInfiniteBus/GenInfiniteBus.cpp
+++ b/examples/Experimental/GenInfiniteBus/GenInfiniteBus.cpp
@@ -1,6 +1,6 @@
-
 #include <iomanip>
 #include <iostream>
+#include <memory>
 
 #include <IpIpoptApplication.hpp>
 #include <IpSolveStatistics.hpp>
@@ -20,49 +20,49 @@ int main()
   using namespace GridKit::Testing;
 
   // Create an infinite bus
-  BaseBus<double, size_t>* bus = new BusSlack<double, size_t>(1.0, 0.0);
+  BusSlack<double, size_t> bus(1.0, 0.0);
 
   // Attach a generator to that bus
-  Generator4<double, size_t>* gen = new Generator4<double, size_t>(bus);
+  Generator4<double, size_t> gen(&bus);
 
   // Create a system model
-  SystemModel<double, size_t>* model = new SystemModel<double, size_t>();
-  model->addBus(bus);
-  model->addComponent(gen);
+  SystemModel<double, size_t> model;
+  model.addBus(&bus);
+  model.addComponent(&gen);
 
   // allocate model components
-  model->allocate();
+  model.allocate();
 
   // Create numerical integrator and configure it for the generator model
-  Ida<double, size_t>* idas = new Ida<double, size_t>(model);
+  Ida<double, size_t> idas(&model);
 
   double t_init  = 0.0;
   double t_final = 15.0;
 
   // setup simulation
-  idas->configureSimulation();
-  idas->configureAdjoint();
-  idas->getDefaultInitialCondition();
-  idas->initializeSimulation(t_init);
-  idas->configureQuadrature();
-  idas->initializeQuadrature();
+  idas.configureSimulation();
+  idas.configureAdjoint();
+  idas.getDefaultInitialCondition();
+  idas.initializeSimulation(t_init);
+  idas.configureQuadrature();
+  idas.initializeQuadrature();
 
   double t_fault = 0.1;
   double t_clear = 0.1;
-  idas->runSimulation(t_fault);
-  idas->saveInitialCondition();
+  idas.runSimulation(t_fault);
+  idas.saveInitialCondition();
   // create initial condition after a fault
   {
-    idas->getSavedInitialCondition();
-    idas->initializeSimulation(t_init);
-    gen->V() = 0.0;
-    idas->runSimulation(t_clear, 20);
-    gen->V() = 1.0;
-    idas->saveInitialCondition();
+    idas.getSavedInitialCondition();
+    idas.initializeSimulation(t_init);
+    gen.V() = 0.0;
+    idas.runSimulation(t_clear, 20);
+    gen.V() = 1.0;
+    idas.saveInitialCondition();
   }
 
   // Set integration time for dynamic constrained optimization
-  idas->setIntegrationTime(t_init, t_final, 100);
+  idas.setIntegrationTime(t_init, t_final, 100);
 
   // Guess initial values of optimization parameters
   double Pm = 1.0;
@@ -90,11 +90,11 @@ int main()
 
   // Create dynamic objective interface to Ipopt solver
   Ipopt::SmartPtr<Ipopt::TNLP> ipoptDynamicObjectiveInterface =
-      new IpoptInterface::DynamicObjective<double, size_t>(idas);
+      new IpoptInterface::DynamicObjective<double, size_t>(&idas);
 
   // Initialize the problem
-  model->param()[0] = Pm;
-  model->param()[1] = Ef;
+  model.param()[0] = Pm;
+  model.param()[1] = Ef;
 
   // Solve the problem
   status = ipoptApp->OptimizeTNLP(ipoptDynamicObjectiveInterface);
@@ -105,26 +105,26 @@ int main()
     // Print result
     std::cout << "\nSucess:\n The problem solved in "
               << ipoptApp->Statistics()->IterationCount() << " iterations!\n"
-              << " Optimal value of Pm = " << model->param()[0] << "\n"
-              << " Optimal value of Ef = " << model->param()[1] << "\n"
+              << " Optimal value of Pm = " << model.param()[0] << "\n"
+              << " Optimal value of Ef = " << model.param()[1] << "\n"
               << " The final value of the objective function G(Pm,Ef) = "
               << ipoptApp->Statistics()->FinalObjective() << "\n\n";
   }
 
   // Create dynamic constraint interface to Ipopt solver
   Ipopt::SmartPtr<Ipopt::TNLP> ipoptDynamicConstraintInterface =
-      new IpoptInterface::DynamicConstraint<double, size_t>(idas);
+      new IpoptInterface::DynamicConstraint<double, size_t>(&idas);
 
   // Store dynamic objective optimization results
-  double* results = new double[model->sizeParams()];
-  for (unsigned i = 0; i < model->sizeParams(); ++i)
+  double* results = new double[model.sizeParams()];
+  for (unsigned i = 0; i < model.sizeParams(); ++i)
   {
-    results[i] = model->param()[i];
+    results[i] = model.param()[i];
   }
 
   // Initialize the problem
-  model->param()[0] = Pm;
-  model->param()[1] = Ef;
+  model.param()[0] = Pm;
+  model.param()[1] = Ef;
 
   // Solve the problem
   status = ipoptApp->OptimizeTNLP(ipoptDynamicConstraintInterface);
@@ -135,17 +135,17 @@ int main()
     // Print result
     std::cout << "\nSucess:\n The problem solved in "
               << ipoptApp->Statistics()->IterationCount() << " iterations!\n"
-              << " Optimal value of Pm = " << model->param()[0] << "\n"
-              << " Optimal value of Ef = " << model->param()[1] << "\n"
+              << " Optimal value of Pm = " << model.param()[0] << "\n"
+              << " Optimal value of Ef = " << model.param()[1] << "\n"
               << " The final value of the objective function G(Pm,Ef) = "
               << ipoptApp->Statistics()->FinalObjective() << "\n\n";
   }
 
   // Compare results of the two optimization methods
   int retval = 0;
-  for (unsigned i = 0; i < model->sizeParams(); ++i)
+  for (unsigned i = 0; i < model.sizeParams(); ++i)
   {
-    if (!isEqual(results[i], model->param()[i], 100 * tol))
+    if (!isEqual(results[i], model.param()[i], 100 * tol))
       --retval;
   }
 
@@ -155,7 +155,5 @@ int main()
   }
 
   delete[] results;
-  delete idas;
-  delete model;
   return 0;
 }
diff --git a/examples/Experimental/ParameterEstimation/ParameterEstimation.cpp b/examples/Experimental/ParameterEstimation/ParameterEstimation.cpp
index 0899dfa4..99598188 100644
--- a/examples/Experimental/ParameterEstimation/ParameterEstimation.cpp
+++ b/examples/Experimental/ParameterEstimation/ParameterEstimation.cpp
@@ -1,6 +1,6 @@
-
 #include <iomanip>
 #include <iostream>
+#include <memory>
 #include <string>
 
 #include "lookup_table.hpp"
@@ -23,58 +23,58 @@ int main()
   using namespace GridKit::Testing;
 
   // Create an infinite bus
-  BaseBus<double, size_t>* bus = new BusSlack<double, size_t>(1.0, 0.0);
+  BusSlack<double, size_t> bus(1.0, 0.0);
 
   // Attach a generator to that bus
-  Generator4Param<double, size_t>* gen = new Generator4Param<double, size_t>(bus);
+  Generator4Param<double, size_t> gen(&bus);
 
   // Create a system model
-  SystemModel<double, size_t>* model = new SystemModel<double, size_t>();
-  model->addBus(bus);
-  model->addComponent(gen);
+  SystemModel<double, size_t> model;
+  model.addBus(&bus);
+  model.addComponent(&gen);
 
   // allocate model components
-  model->allocate();
+  model.allocate();
 
   // Create numerical integrator and configure it for the generator model
-  Ida<double, size_t>* idas = new Ida<double, size_t>(model);
+  Ida<double, size_t> idas(&model);
 
   double t_init  = -1.0;
   double t_final = -1.0;
 
   std::istringstream input_data(lookup_table);
-  GridKit::setLookupTable(gen->getLookupTable(), input_data, t_init, t_final);
+  GridKit::setLookupTable(gen.getLookupTable(), input_data, t_init, t_final);
 
   std::cout << "Performing parameter estimation with respect to data\nfrom "
             << "t_init = " << t_init << " to t_final = " << t_final << "\n";
 
   // setup simulation
-  idas->configureSimulation();
-  idas->configureAdjoint();
-  idas->getDefaultInitialCondition();
-  idas->initializeSimulation(t_init);
-  idas->configureQuadrature();
-  idas->initializeQuadrature();
+  idas.configureSimulation();
+  idas.configureAdjoint();
+  idas.getDefaultInitialCondition();
+  idas.initializeSimulation(t_init);
+  idas.configureQuadrature();
+  idas.initializeQuadrature();
 
   double t_fault = 0.1;
   double t_clear = 0.1;
-  idas->runSimulation(t_fault);
-  idas->saveInitialCondition();
+  idas.runSimulation(t_fault);
+  idas.saveInitialCondition();
   // create initial condition after a fault
   {
-    idas->getSavedInitialCondition();
-    idas->initializeSimulation(t_init);
-    gen->V() = 0.0;
-    idas->runSimulation(t_clear, 20);
-    gen->V() = 1.0;
-    idas->saveInitialCondition();
+    idas.getSavedInitialCondition();
+    idas.initializeSimulation(t_init);
+    gen.V() = 0.0;
+    idas.runSimulation(t_clear, 20);
+    gen.V() = 1.0;
+    idas.saveInitialCondition();
   }
 
   // Set integration time for dynamic constrained optimization
-  idas->setIntegrationTime(t_init, t_final, 100);
+  idas.setIntegrationTime(t_init, t_final, 100);
 
   // Guess value of inertia coefficient
-  model->param()[0] = 3.0;
+  model.param()[0] = 3.0;
 
   // Create an instance of the IpoptApplication
   Ipopt::SmartPtr<Ipopt::IpoptApplication> ipoptApp = IpoptApplicationFactory();
@@ -98,7 +98,7 @@ int main()
 
   // Create dynamic objective interface to Ipopt solver
   Ipopt::SmartPtr<Ipopt::TNLP> ipoptDynamicObjectiveInterface =
-      new IpoptInterface::DynamicObjective<double, size_t>(idas);
+      new IpoptInterface::DynamicObjective<double, size_t>(&idas);
 
   // Solve the problem
   status = ipoptApp->OptimizeTNLP(ipoptDynamicObjectiveInterface);
@@ -109,24 +109,24 @@ int main()
     // Print result
     std::cout << "\nSucess:\n The problem solved in "
               << ipoptApp->Statistics()->IterationCount() << " iterations!\n"
-              << " Optimal value of H = " << model->param()[0] << "\n"
+              << " Optimal value of H = " << model.param()[0] << "\n"
               << " The final value of the objective function G(H) = "
               << ipoptApp->Statistics()->FinalObjective() << "\n\n";
   }
 
   // Store dynamic objective optimization results
-  double* results = new double[model->sizeParams()];
-  for (unsigned i = 0; i < model->sizeParams(); ++i)
+  double* results = new double[model.sizeParams()];
+  for (unsigned i = 0; i < model.sizeParams(); ++i)
   {
-    results[i] = model->param()[i];
+    results[i] = model.param()[i];
   }
 
   // Guess value of inertia coefficient
-  model->param()[0] = 3.0;
+  model.param()[0] = 3.0;
 
   // Create dynamic constraint interface to Ipopt solver
   Ipopt::SmartPtr<Ipopt::TNLP> ipoptDynamicConstraintInterface =
-      new IpoptInterface::DynamicConstraint<double, size_t>(idas);
+      new IpoptInterface::DynamicConstraint<double, size_t>(&idas);
 
   // Solve the problem
   status = ipoptApp->OptimizeTNLP(ipoptDynamicConstraintInterface);
@@ -137,16 +137,16 @@ int main()
     // Print result
     std::cout << "\nSucess:\n The problem solved in "
               << ipoptApp->Statistics()->IterationCount() << " iterations!\n"
-              << " Optimal value of H = " << model->param()[0] << "\n"
+              << " Optimal value of H = " << model.param()[0] << "\n"
               << " The final value of the objective function G(H) = "
               << ipoptApp->Statistics()->FinalObjective() << "\n\n";
   }
 
   // Compare results of the two optimization methods
   int retval = 0;
-  for (unsigned i = 0; i < model->sizeParams(); ++i)
+  for (unsigned i = 0; i < model.sizeParams(); ++i)
   {
-    if (!isEqual(results[i], model->param()[i], 100 * tol))
+    if (!isEqual(results[i], model.param()[i], 100 * tol))
       --retval;
   }
 
@@ -156,7 +156,5 @@ int main()
   }
 
   delete[] results;
-  delete idas;
-  delete model;
   return retval;
 }
diff --git a/src/Solver/Dynamic/Ida.cpp b/src/Solver/Dynamic/Ida.cpp
index 01035a9d..bc71866f 100644
--- a/src/Solver/Dynamic/Ida.cpp
+++ b/src/Solver/Dynamic/Ida.cpp
@@ -425,15 +425,23 @@ namespace AnalysisManager
       retval = IDASetMaxNumStepsB(solver_, backwardID_, 2000);
       checkOutput(retval, "IDASetMaxNumSteps");
 
-      // Set up linear solver
-      JacobianMatB_ = SUNDenseMatrix(static_cast<sunindextype>(model_->size()),
-                                     static_cast<sunindextype>(model_->size()),
-                                     context_);
-      checkAllocation((void*) JacobianMatB_, "SUNDenseMatrix");
+      // Allocate Jacobian matrix, if not already
+      if (JacobianMatB_ == nullptr)
+      {
+        JacobianMatB_ = SUNDenseMatrix(static_cast<sunindextype>(model_->size()),
+                                       static_cast<sunindextype>(model_->size()),
+                                       context_);
+        checkAllocation((void*) JacobianMatB_, "SUNDenseMatrix");
+      }
 
-      linearSolverB_ = SUNLinSol_Dense(yyB_, JacobianMatB_, context_);
-      checkAllocation((void*) linearSolverB_, "SUNLinSol_Dense");
+      // Allocate linear solver, if not already
+      if (linearSolverB_ == nullptr)
+      {
+        linearSolverB_ = SUNLinSol_Dense(yyB_, JacobianMatB_, context_);
+        checkAllocation((void*) linearSolverB_, "SUNLinSol_Dense");
+      }
 
+      // Setup linear solver (only dense supported at this time)
       retval = IDASetLinearSolverB(solver_, backwardID_, linearSolverB_, JacobianMatB_);
       checkOutput(retval, "IDASetLinearSolverB");