diff --git a/CHANGELOG.md b/CHANGELOG.md
index 04aa4071..ca382a44 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -40,6 +40,7 @@
 - Added `GridKitDocs` target for Doxygen documentation.
 - Added a header file defining common math functions (e.g., sigmoid) to be used throughout the code.
 - Added capability to print monitored variables in multiple formats, triggered from `Ida::runSimulation`.
+- Added IDA statistics object which can be accumulated over multiple simulations.
 
 ## v0.1
 
diff --git a/GridKit/Solver/Dynamic/Ida.cpp b/GridKit/Solver/Dynamic/Ida.cpp
index baa303a5..a05ba1ca 100644
--- a/GridKit/Solver/Dynamic/Ida.cpp
+++ b/GridKit/Solver/Dynamic/Ida.cpp
@@ -3,6 +3,7 @@
 
 #include <iomanip>
 #include <iostream>
+#include <sstream>
 
 #include <idas/idas.h>
 #include <idas/idas_ls.h>
@@ -968,6 +969,76 @@ namespace AnalysisManager
       checkOutput(retval, "IDAPrintAllStats");
     }
 
+    /**
+     * @brief Accumulate another stats object into this one, allowing for stats to be kept
+     *        across multiple simulations with IDA.
+     */
+    IdaStats& IdaStats::operator+=(const IdaStats& other)
+    {
+      num_steps_                       += other.num_steps_;
+      num_residual_evals_              += other.num_residual_evals_;
+      num_linear_decompositions_       += other.num_linear_decompositions_;
+      num_error_test_fails_            += other.num_error_test_fails_;
+      num_nonlinear_iters_             += other.num_nonlinear_iters_;
+      num_nonlinear_convergence_fails_ += other.num_nonlinear_convergence_fails_;
+
+      return *this;
+    }
+
+    /**
+     * @brief Generate a string containing all of the stats in a formatted report.
+     *        All columns are aligned. To change the width of a column,
+     *        modify `label_width` or `stat_width`.
+     */
+    std::string IdaStats::report() const
+    {
+      int               label_width = 30;
+      int               stat_width  = 12;
+      std::stringstream out;
+
+      out << std::setw(label_width) << "Steps" << " : " << std::setw(stat_width) << num_residual_evals_ << '\n'
+          << std::setw(label_width) << "Residual evals" << " : " << std::setw(stat_width) << num_linear_decompositions_ << '\n'
+          << std::setw(label_width) << "Linear decompositions" << " : " << std::setw(stat_width) << num_linear_decompositions_ << '\n'
+          << std::setw(label_width) << "Error test failures" << " : " << std::setw(stat_width) << num_error_test_fails_ << '\n'
+          << std::setw(label_width) << "Nonlinear iterations" << " : " << std::setw(stat_width) << num_nonlinear_iters_ << '\n'
+          << std::setw(label_width) << "Nonlinear convergence failures" << " : " << std::setw(stat_width) << num_nonlinear_convergence_fails_;
+
+      return out.str();
+    }
+
+    /**
+     * @brief Construct and return an `IdaStats` object containing the statistics of the current IDA workspace.
+     *        Several statistics returned by IDA are ignored because they are about the current state of IDA,
+     *        rather than about the simulation at large.
+     */
+    template <class ScalarT, typename IdxT>
+    IdaStats Ida<ScalarT, IdxT>::getStats() const
+    {
+      IdaStats stats;
+
+      // Dummies for ignoring stats
+      int         dummy;
+      sunrealtype dummy2;
+
+      int retval = IDAGetIntegratorStats(solver_,
+                                         &stats.num_steps_,
+                                         &stats.num_residual_evals_,
+                                         &stats.num_linear_decompositions_,
+                                         &stats.num_error_test_fails_,
+                                         &dummy,
+                                         &dummy,
+                                         &dummy2,
+                                         &dummy2,
+                                         &dummy2,
+                                         &dummy2);
+      checkOutput(retval, "IDAGetIntegratorStats");
+
+      retval = IDAGetNonlinSolvStats(solver_, &stats.num_nonlinear_iters_, &stats.num_nonlinear_convergence_fails_);
+      checkOutput(retval, "IDAGetNonlinSolvStats");
+
+      return stats;
+    }
+
     /**
      * @brief Check SUNDIALS allocation
      *
diff --git a/GridKit/Solver/Dynamic/Ida.hpp b/GridKit/Solver/Dynamic/Ida.hpp
index 8e1715e5..650e2885 100644
--- a/GridKit/Solver/Dynamic/Ida.hpp
+++ b/GridKit/Solver/Dynamic/Ida.hpp
@@ -24,6 +24,19 @@ namespace AnalysisManager
 {
   namespace Sundials
   {
+    struct IdaStats
+    {
+      long int num_steps_                       = 0;
+      long int num_residual_evals_              = 0;
+      long int num_linear_decompositions_       = 0;
+      long int num_error_test_fails_            = 0;
+      long int num_nonlinear_iters_             = 0;
+      long int num_nonlinear_convergence_fails_ = 0;
+
+      IdaStats&   operator+=(const IdaStats& other);
+      std::string report() const;
+    };
+
     template <class ScalarT, typename IdxT>
     class Ida : public DynamicSolver<ScalarT, IdxT>
     {
@@ -115,6 +128,8 @@ namespace AnalysisManager
       void printSpecial(RealT t, N_Vector x);
       void printFinalStats();
 
+      IdaStats getStats() const;
+
     private:
       static int Residual(RealT    t,
                           N_Vector yy,
@@ -188,8 +203,8 @@ namespace AnalysisManager
       static void copyVec(const std::vector<bool>& x, N_Vector y);
 
       // int check_flag(void *flagvalue, const char *funcname, int opt);
-      inline void checkAllocation(void* v, const char* functionName);
-      inline void checkOutput(int retval, const char* functionName);
+      static void checkAllocation(void* v, const char* functionName);
+      static void checkOutput(int retval, const char* functionName);
     };
 
     /// Simple exception to use within Ida class.
diff --git a/examples/PowerElectronics/RLCircuit/RLCircuit.cpp b/examples/PowerElectronics/RLCircuit/RLCircuit.cpp
index dce13668..44804aa7 100644
--- a/examples/PowerElectronics/RLCircuit/RLCircuit.cpp
+++ b/examples/PowerElectronics/RLCircuit/RLCircuit.cpp
@@ -134,5 +134,7 @@ int main(int /* argc */, char const** /* argv */)
     std::cout << abs((yfinial[i] - yexact[i]) / yexact[i]) << "\n";
   }
 
+  std::cerr << idas.getStats().report() << '\n';
+
   return 0;
 }