diff --git a/libnufi_lib.a b/libnufi_lib.a
index f66070d..0b9082d 100644
Binary files a/libnufi_lib.a and b/libnufi_lib.a differ
diff --git a/nufi/fields.h b/nufi/fields.h
index 9268c47..84ed5a7 100644
--- a/nufi/fields.h
+++ b/nufi/fields.h
@@ -9,52 +9,54 @@
 
 using namespace dealii;
 
-inline std::vector<int> Indices_of_points(const std::vector<double> &points, double x_min, double x_max, double dx, int grid_type=0)
-{ 
-    // grid type:
-    // 0 => uniform
-    // 1 => non uniform (TODO)
-    
-    if (dx <= 0.0) {
-        throw std::invalid_argument("dx must be positive");
-    }
-    if (x_max <= x_min) {
-        throw std::invalid_argument("x_max must be > x_min");
-    }
-
-    std::vector<int> indices;
-    indices.reserve(points.size());
-
-    switch (grid_type) {
-        case 0:
-            {
-            const double L = x_max - x_min;
-            const int N = std::floor(L/dx);
-
-
-            for (double x : points) //GPT loop, to check
-            {
-                x-= x_min;
-                x = x - L * std::floor(x/L);
-
-                int i = static_cast<int>(std::floor(x / dx));
-
-                // safety: handle rare edge case due to floating precision
-                if (i == N) i = 0;
-
-                indices.push_back(i);
-            }
-            }
-        case 1:
-            {
-            throw std::invalid_argument("Case for non uniform grid is not completed");
-            }
-        default:
-            throw std::invalid_argument("Invalid grid_type argument");
-
-    }
-    return indices;
-}
+// inline std::vector<int> Indices_of_points(const std::vector<double> &points,
+// double x_min, double x_max, double dx, int grid_type=0)
+// {
+//     // grid type:
+//     // 0 => uniform
+//     // 1 => non uniform (TODO)
+//
+//     if (dx <= 0.0) {
+//         throw std::invalid_argument("dx must be positive");
+//     }
+//     if (x_max <= x_min) {
+//         throw std::invalid_argument("x_max must be > x_min");
+//     }
+//
+//     std::vector<int> indices;
+//     indices.reserve(points.size());
+//
+//     switch (grid_type) {
+//         case 0:
+//             {
+//             const double L = x_max - x_min;
+//             const int N = std::floor(L/dx);
+//
+//
+//             for (double x : points) //GPT loop, to check
+//             {
+//                 x-= x_min;
+//                 x = x - L * std::floor(x/L);
+//
+//                 int i = static_cast<int>(std::floor(x / dx));
+//
+//                 // safety: handle rare edge case due to floating precision
+//                 if (i == N) i = 0;
+//
+//                 indices.push_back(i);
+//             }
+//             }
+//         case 1:
+//             {
+//             throw std::invalid_argument("Case for non uniform grid is not
+//             completed");
+//             }
+//         default:
+//             throw std::invalid_argument("Invalid grid_type argument");
+//
+//     }
+//     return indices;
+// }
 
 inline double f0(const double x, const double v,
                  const double eps = Parameters::EPS,
@@ -81,11 +83,18 @@ inline double compute_rho(const double x,
   return 1.0 - integral;
 }
 
-double eval(double x, const PoissonProblem<1> &poisson) noexcept {
-    return poisson.evaluate_potential(Point<1>(x));
+inline double eval(double x, const PoissonProblem<1> &poisson,
+                   const Vector<double> &solution) noexcept {
+
+  x -= Parameters::X_DOMAIN_LEFT;
+
+  x = x - Parameters::LX * std::floor(x * Parameters::LX_INV);
+  return eval_point<1>(poisson.get_mapping(), poisson.get_dof_handler(),
+                       solution, Point<1>(x));
 }
 
 inline double integral_space_vector(const PoissonProblem<1> &poisson,
+                                    const Vector<double> &solution,
                                     double dx = Parameters::PLOT_DX,
                                     size_t Nx = Parameters::PLOT_NX) {
   double integral = 0.0;
@@ -93,12 +102,13 @@ inline double integral_space_vector(const PoissonProblem<1> &poisson,
 #pragma omp parallel for reduction(+ : integral)
   for (size_t i = 0; i < Nx; ++i) {
     double x = xmin + i * dx;
-    integral += eval(x, poisson);
+    integral += eval(x, poisson, solution);
   }
   return integral * dx;
 };
 
 inline double integral_space_vector_squared(const PoissonProblem<1> &poisson,
+                                            const Vector<double> &solution,
                                             double dx = Parameters::PLOT_DX,
                                             size_t Nx = Parameters::PLOT_NX) {
   double integral = 0.0;
@@ -106,46 +116,46 @@ inline double integral_space_vector_squared(const PoissonProblem<1> &poisson,
 #pragma omp parallel for reduction(+ : integral)
   for (size_t i = 0; i < Nx; ++i) {
     double x = xmin + i * dx;
-    double val = eval(x, poisson);
+    double val = eval(x, poisson, solution);
     integral += val * val;
   }
   return integral * dx;
 };
 
-class Gradient {
-public:
-  Gradient(double xmin, double xmax, unsigned int Nx)
-      : xmin_(xmin), xmax_(xmax), Nx_(Nx) {
-    if (xmax_ <= xmin_) {
-      throw std::invalid_argument("xmax must be greater than xmin");
-    }
-  }
-
-  std::vector<double> compute(const std::vector<double> &values) const {
-    size_t n = values.size();
-    if (n < 2) {
-      throw std::invalid_argument("Need at least 2 points");
-    }
-
-    std::vector<double> grad(n);
-
-    double dx = (xmax_ - xmin_) / (n - 1);
-    // periodic boundaries
-    grad[0] = -(values[1] - values[n - 1]) / (2.0 * dx);
-    grad[n - 1] = -(values[0] - values[n - 2]) / (2.0 * dx);
-
-    for (size_t i = 1; i < n - 1; ++i) {
-      grad[i] = -(values[i + 1] - values[i - 1]) / (2.0 * dx);
-    }
-
-    return grad;
-  }
-
-private:
-  double xmin_;
-  double xmax_;
-  [[maybe_unused]] unsigned int Nx_;
-};
+// class Gradient {
+// public:
+//   Gradient(double xmin, double xmax, unsigned int Nx)
+//       : xmin_(xmin), xmax_(xmax), Nx_(Nx) {
+//     if (xmax_ <= xmin_) {
+//       throw std::invalid_argument("xmax must be greater than xmin");
+//     }
+//   }
+//
+//   std::vector<double> compute(const std::vector<double> &values) const {
+//     size_t n = values.size();
+//     if (n < 2) {
+//       throw std::invalid_argument("Need at least 2 points");
+//     }
+//
+//     std::vector<double> grad(n);
+//
+//     double dx = (xmax_ - xmin_) / (n - 1);
+//     // periodic boundaries
+//     grad[0] = -(values[1] - values[n - 1]) / (2.0 * dx);
+//     grad[n - 1] = -(values[0] - values[n - 2]) / (2.0 * dx);
+//
+//     for (size_t i = 1; i < n - 1; ++i) {
+//       grad[i] = -(values[i + 1] - values[i - 1]) / (2.0 * dx);
+//     }
+//
+//     return grad;
+//   }
+//
+// private:
+//   double xmin_;
+//   double xmax_;
+//   [[maybe_unused]] unsigned int Nx_;
+// };
 
 template <int dim>
 class ChargeDensity : public Function<dim> // only uses f0
diff --git a/nufi/nufi_solver.h b/nufi/nufi_solver.h
index 522b0b5..9af3813 100644
--- a/nufi/nufi_solver.h
+++ b/nufi/nufi_solver.h
@@ -2,32 +2,35 @@
 #define NUFI_SOLVER_H
 
 #include <boost/qvm/mat_access.hpp>
-#include <vector>
 #include <cmath>
 #include <deal.II/base/point.h>
 #include <deal.II/base/tensor.h>
 #include <deal.II/numerics/fe_field_function.h>
+#include <vector>
 
+#include "nufi/fields.h" //dont remove
 #include "nufi/parameters.h"
 #include "nufi/poisson_problem.h"
-#include "nufi/fields.h" //dont remove
 
 using namespace dealii;
 
-class NuFISolver
-{
+class NuFISolver {
 public:
   NuFISolver();
 
   void run();
-  double eval_rho(unsigned int n, double x, const PoissonProblem<1> &poisson, unsigned int Nv = Parameters::NV) const;
-  double eval_ftilda(unsigned int n, double x, double u, const PoissonProblem<1> &poisson) const;
-  double eval_f(unsigned int n, double x, double u, const PoissonProblem<1> &poisson) const;
+  double eval_rho(unsigned int n, double x, const PoissonProblem<1> &poisson,
+                  const std::vector<Vector<double>> &phi_history,
+                  unsigned int Nv = Parameters::NV) const;
+  double eval_ftilda(unsigned int n, double x, double u,
+                     const PoissonProblem<1> &poisson,
+                     const std::vector<Vector<double>> &phi_history) const;
+  double eval_f(unsigned int n, double x, double u,
+                const PoissonProblem<1> &poisson,
+                const std::vector<Vector<double>> &phi_history) const;
 
 private:
-
-
-  unsigned int Nt = std::floor(Parameters::TMAX/Parameters::DT);
+  unsigned int Nt = std::floor(Parameters::TMAX / Parameters::DT);
   unsigned int Nx = Parameters::CALC_NX;
 
   double Lx = Parameters::LX;
@@ -40,36 +43,33 @@ private:
   unsigned int order;
 
   PoissonProblem<1> poisson;
-
 };
 
-template<unsigned int dim>
-class ChargeDensity_NuFI : public Function<dim>
-{
-  public:
-    ChargeDensity_NuFI(const double *rho_values, unsigned int Nx)
+template <unsigned int dim> class ChargeDensity_NuFI : public Function<dim> {
+public:
+  ChargeDensity_NuFI(const double *rho_values, unsigned int Nx)
       : Function<dim>(), rho(rho_values), Nx(Nx) {}
 
-    virtual double value(const Point<dim> &p,
-                         [[maybe_unused]] const unsigned int component = 0) const override
-    {
-      const double x = p[0];
+  virtual double
+  value(const Point<dim> &p,
+        [[maybe_unused]] const unsigned int component = 0) const override {
+    const double x = p[0];
 
-      // Map x -> grid index
-      const double L = Parameters::LX;
-      const double dx = L / (Nx-1);
+    // Map x -> grid index
+    const double L = Parameters::LX;
+    const double dx = L / (Nx - 1);
 
-      int i = static_cast<int>(std::floor((x - Parameters::X_DOMAIN_LEFT) / dx));
+    int i = static_cast<int>(std::floor((x - Parameters::X_DOMAIN_LEFT) / dx));
 
-      // periodic wrap
-      i = (i % Nx + Nx) % Nx;
+    // periodic wrap
+    i = (i % Nx + Nx) % Nx;
 
-      return rho[i];
-    }
+    return rho[i];
+  }
 
-  private:
-    const double *rho;
-    const unsigned int Nx;
+private:
+  const double *rho;
+  const unsigned int Nx;
 };
 
 #endif
diff --git a/nufi/parameters.h b/nufi/parameters.h
index 0049d18..a457124 100644
--- a/nufi/parameters.h
+++ b/nufi/parameters.h
@@ -4,42 +4,41 @@
 #include <cmath>
 #include <cstdlib>
 
-namespace Parameters
-{
-  constexpr unsigned int DIMENSION = 1;
+namespace Parameters {
+constexpr unsigned int DIMENSION = 1;
 
-  constexpr double X_DOMAIN_LEFT  = 0.0;
-  constexpr double X_DOMAIN_RIGHT = 4*M_PI;
-  constexpr double LX = std::abs(X_DOMAIN_RIGHT- X_DOMAIN_LEFT);
-  constexpr double LX_INV = 1/LX;
+constexpr double X_DOMAIN_LEFT = 0.0;
+constexpr double X_DOMAIN_RIGHT = 4 * M_PI;
+constexpr double LX = std::abs(X_DOMAIN_RIGHT - X_DOMAIN_LEFT);
+constexpr double LX_INV = 1 / LX;
 
-  constexpr size_t CALC_NX = 256;
-  constexpr double CALC_DX = LX/CALC_NX;
+constexpr size_t CALC_NX = 128;
+constexpr double CALC_DX = LX / CALC_NX;
 
-  constexpr double V_DOMAIN_LEFT = -10.;
-  constexpr double V_DOMAIN_RIGHT = 10.;
+constexpr double V_DOMAIN_LEFT = -10.;
+constexpr double V_DOMAIN_RIGHT = 10.;
 
-  constexpr unsigned int NV = 256;
-  constexpr double DV = std::abs(V_DOMAIN_RIGHT - V_DOMAIN_LEFT)/NV;
+constexpr unsigned int NV = 128;
+constexpr double DV = std::abs(V_DOMAIN_RIGHT - V_DOMAIN_LEFT) / NV;
 
-  // deal.ii options
-  constexpr unsigned int GLOBAL_REFINEMENT = 8;
-  constexpr unsigned int FE_DEGREE = 3;
-  constexpr unsigned int CONVERGENCE_ITERATIONS = 5000;
-  constexpr double CONVERGENCE_LIMIT = 1e-8;
+// deal.ii options
+constexpr unsigned int GLOBAL_REFINEMENT = 6;
+constexpr unsigned int FE_DEGREE = 2;
+constexpr unsigned int CONVERGENCE_ITERATIONS = 5000;
+constexpr double CONVERGENCE_LIMIT = 1e-8;
 
-  constexpr double EPS = 0.01;
-  constexpr double WAVE_NR = 0.5;
-  constexpr double F0_FACTOR = 0.39894228040143267793994; // 1/sqrt(2pi)
+constexpr double EPS = 0.01;
+constexpr double WAVE_NR = 0.5;
+constexpr double F0_FACTOR = 0.39894228040143267793994; // 1/sqrt(2pi)
 
-  // NUFI options
-  constexpr double DT=1./4.;
-  constexpr unsigned int TMAX = 50;
+// NUFI options
+constexpr double DT = 1. / 10.;
+constexpr unsigned int TMAX = 100;
 
-  //Plotting options
-  constexpr int PLOT_FREQUENCY = 20;
-  constexpr size_t PLOT_NX = 256;
-  constexpr double PLOT_DX = LX/PLOT_NX;
-}
+// Plotting options
+constexpr int PLOT_FREQUENCY = 10;
+constexpr size_t PLOT_NX = CALC_NX;
+constexpr double PLOT_DX = LX / PLOT_NX;
+} // namespace Parameters
 
 #endif
diff --git a/nufi/poisson_problem.h b/nufi/poisson_problem.h
index 3de0c54..cae16a3 100644
--- a/nufi/poisson_problem.h
+++ b/nufi/poisson_problem.h
@@ -12,6 +12,7 @@
 #include <deal.II/base/tensor.h>
 #include <deal.II/base/utilities.h>
 
+#include <deal.II/fe/mapping_q.h>
 #include <deal.II/lac/affine_constraints.h>
 #include <deal.II/lac/dynamic_sparsity_pattern.h>
 #include <deal.II/lac/full_matrix.h>
@@ -61,17 +62,13 @@ public:
 
   const Vector<double> &get_solution() const { return solution; }
   const DoFHandler<dim> &get_dof_handler() const { return dof_handler; }
+  const MappingQ<dim> &get_mapping() const { return mapping; }
 
   std::vector<double> sample_electric_field(double x_min, double x_max,
                                             unsigned int Nx);
   std::vector<double> sample_electric_potential(double x_min, double x_max,
                                                 unsigned int Nx);
 
-  double evaluate_potential(const Point<dim> &p) const
-{
-  return fe_field_function->value(p);
-}
-
 private:
   void create_mesh();
   void setup_system();
@@ -94,7 +91,7 @@ private:
 
   MappingQ<dim> mapping;
 
-  std::unique_ptr<Functions::FEFieldFunction<dim>> fe_field_function;
+  // std::unique_ptr<Functions::FEFieldFunction<dim>> fe_field_function;
 };
 
 // Utilities
@@ -177,8 +174,8 @@ PoissonProblem<dim>::sample_electric_potential(double x_min, double x_max,
 // template <int dim> std::vector<double> eval_solution_on_points(
 //         const std::vector<Vector<double>> &solutions,
 //         const unsigned int n,
-//         const std::vector<Point<dim>> &points, // need to be in [x_min, x_max]. I think....
-//         const std::vector<unsigned int> &cell_indices,
+//         const std::vector<Point<dim>> &points, // need to be in [x_min,
+//         x_max]. I think.... const std::vector<unsigned int> &cell_indices,
 //         const DoFHandler<dim> &dof_handler,
 //         const MappingQ<dim> &mapping)
 // {
@@ -219,7 +216,8 @@ PoissonProblem<dim>::sample_electric_potential(double x_min, double x_max,
 //         for (unsigned int id : point_ids)
 //             cell_points.push_back(points[id]);
 //
-//         std::vector<types::global_dof_index> indices(dof_handler.get_fe().n_dofs_per_cell());
+//         std::vector<types::global_dof_index>
+//         indices(dof_handler.get_fe().n_dofs_per_cell());
 //         cell->get_dof_indices(indices);
 //
 //         for (unsigned int i=0;i<indices.size();++i)
@@ -240,13 +238,8 @@ PoissonProblem<dim>::sample_electric_potential(double x_min, double x_max,
 template <int dim>
 double eval_point(const Mapping<dim> &mapping,
                   const DoFHandler<dim> &dof_handler,
-                  const Vector<double> &solution,
-                  const Point<dim> &point)
-{
-    return VectorTools::point_value<dim>(mapping,
-                                         dof_handler,
-                                         solution,
-                                         point);
+                  const Vector<double> &solution, const Point<dim> &point) {
+  return VectorTools::point_value<dim>(mapping, dof_handler, solution, point);
 }
 
 // dealii Poisson
@@ -288,7 +281,7 @@ template <int dim> void PoissonProblem<dim>::setup_system() {
       gauge_dof = i;
       break;
     }
-
+  }
 
   Assert(gauge_dof != numbers::invalid_dof_index,
          ExcMessage("No unconstrained DoF found for gauge fixing."));
@@ -307,11 +300,9 @@ template <int dim> void PoissonProblem<dim>::setup_system() {
   solution.reinit(dof_handler.n_dofs());
   system_rhs.reinit(dof_handler.n_dofs());
 
-  fe_field_function =
-  std::make_unique<Functions::FEFieldFunction<dim>>(
-    dof_handler, solution, mapping);
-  }
-
+  // fe_field_function =
+  // std::make_unique<Functions::FEFieldFunction<dim>>(
+  //   dof_handler, solution, mapping);
 }
 /* (Mine)
 template <int dim>
@@ -393,6 +384,8 @@ void PoissonProblem<dim>::assemble_system()
 // Paul's, mine's above
 
 template <int dim> void PoissonProblem<dim>::assemble_system() {
+  Assert(system_matrix.m() == dof_handler.n_dofs(),
+         ExcMessage("Matrix not initialized correctly"));
   system_matrix = 0;
   system_rhs = 0;
 
@@ -446,11 +439,6 @@ template <int dim> void PoissonProblem<dim>::solve() {
   // solver.solve(system_matrix, solution, system_rhs, preconditioner);
   solver.solve(system_matrix, solution, system_rhs, PreconditionIdentity());
   constraints.distribute(solution);
-
-
-  fe_field_function =
-  std::make_unique<Functions::FEFieldFunction<dim>>(
-    dof_handler, solution, mapping);
 }
 
 template <int dim> void PoissonProblem<dim>::initialize() {
diff --git a/nufi/save_results.h b/nufi/save_results.h
index 981158c..7e47ae7 100644
--- a/nufi/save_results.h
+++ b/nufi/save_results.h
@@ -1,29 +1,26 @@
 #ifndef SAVE_RESULTS_H
 #define SAVE_RESULTS_H
 
-#include <string>
 #include "nufi/nufi_solver.h"
 #include "nufi/poisson_problem.h"
+#include <deal.II/lac/vector.h>
+#include <string>
 
+void save_f(const NuFISolver &solver, unsigned int n,
+            const PoissonProblem<1> &poisson,
+            const std::vector<Vector<double>> &phi_history, unsigned int Nx_out,
+            unsigned int Nv_out, const std::string &filename);
 
-void save_f( const NuFISolver &solver,
-    unsigned int n,
-                                    const PoissonProblem<1> &poisson,
-                                    unsigned int Nx_out,
-                                    unsigned int Nv_out,
-                                    const std::string &filename);
+void save_rho(const NuFISolver &solver, unsigned int n,
+              const PoissonProblem<1> &poisson,
+              const std::vector<Vector<double>> &phi_history,
+              unsigned int Nx_out, const std::string &filename);
 
-void save_rho(const NuFISolver &solver,
-    unsigned int n,
-                                    const PoissonProblem<1> &poisson,
-                                    unsigned int Nx_out,
-                                    const std::string &filename);
+void save_Efield(unsigned int n, const PoissonProblem<1> &poisson,
+                 const std::vector<Vector<double>> &phi_history,
+                 unsigned int Nx_out, const std::string &filename);
 
-void save_Efield(unsigned int n,
-                                    const PoissonProblem<1> &poisson,
-                                    unsigned int Nx_out,
-                                    const std::string &filename);
+void save_space_vector(const std::vector<double> &vals,
+                       const std::string &filename, size_t it);
 
-void save_space_vector(const std::vector<double>& vals, const std::string& filename, size_t it);
-
-#endif 
+#endif
diff --git a/nufi/stopwatch.h b/nufi/stopwatch.h
index 1d2b5ca..2c6bbe8 100644
--- a/nufi/stopwatch.h
+++ b/nufi/stopwatch.h
@@ -3,35 +3,27 @@
 
 #include <chrono>
 
-template <typename real>
-class stopwatch
-{
+template <typename real> class stopwatch {
 public:
-	void reset();
-	real elapsed();
+  void reset();
+  real elapsed();
 
 private:
-	using clock = std::chrono::high_resolution_clock;
-	clock::time_point t0 { clock::now() };
+  using clock = std::chrono::high_resolution_clock;
+  clock::time_point t0{clock::now()};
 };
 
-
-template <typename real> inline
-void stopwatch<real>::reset()
-{
-	t0 = clock::now();
+template <typename real> inline void stopwatch<real>::reset() {
+  t0 = clock::now();
 }
 
-template <typename real> inline
-real stopwatch<real>::elapsed()
-{
-	using seconds = std::chrono::duration<real,std::ratio<1,1>>;
+template <typename real> inline real stopwatch<real>::elapsed() {
+  using seconds = std::chrono::duration<real, std::ratio<1, 1>>;
 
-	auto tnow = clock::now();
-	auto duration = std::chrono::duration_cast<seconds>( tnow - t0 );
+  auto tnow = clock::now();
+  auto duration = std::chrono::duration_cast<seconds>(tnow - t0);
 
-	return duration.count();
+  return duration.count();
 }
 
 #endif // STOPWATCH_H
-
diff --git a/src/main.cc b/src/main.cc
index 7e18060..66caca3 100644
--- a/src/main.cc
+++ b/src/main.cc
@@ -1,44 +1,35 @@
-#include <iostream>
 #include <filesystem>
+#include <iostream>
 
 #include <nufi/nufi_solver.h>
 
-void clear_results_directory(const std::string &dir)
-{
-    if (!std::filesystem::exists(dir) || !std::filesystem::is_directory(dir))
-        return;
+void clear_results_directory(const std::string &dir) {
+  if (!std::filesystem::exists(dir) || !std::filesystem::is_directory(dir))
+    return;
 
-    for (const auto &entry : std::filesystem::directory_iterator(dir))
-    {
-        if (std::filesystem::is_regular_file(entry))
-        {
-            std::filesystem::remove(entry.path());
-            std::cout << "Deleted: " << entry.path() << '\n';
-        }
+  for (const auto &entry : std::filesystem::directory_iterator(dir)) {
+    if (std::filesystem::is_regular_file(entry)) {
+      std::filesystem::remove(entry.path());
+      std::cout << "Deleted: " << entry.path() << '\n';
     }
+  }
 }
 
-int main()
-{
-  #include <omp.h>
-std::cout << "Threads: " << omp_get_max_threads() << "\n";
-    try
-    {
-        clear_results_directory("results");
+int main() {
+#include <omp.h>
+  std::cout << "Threads: " << omp_get_max_threads() << "\n";
+  try {
+    clear_results_directory("results");
 
-        NuFISolver solver;
-        solver.run();
-    }
-    catch (const std::exception &exc)
-    {
-        std::cerr << "\nException:\n" << exc.what() << "\n";
-        return 1;
-    }
-    catch (...)
-    {
-        std::cerr << "\nUnknown exception!\n";
-        return 1;
-    }
+    NuFISolver solver;
+    solver.run();
+  } catch (const std::exception &exc) {
+    std::cerr << "\nException:\n" << exc.what() << "\n";
+    return 1;
+  } catch (...) {
+    std::cerr << "\nUnknown exception!\n";
+    return 1;
+  }
 
-    return 0;
+  return 0;
 }
diff --git a/src/nufi_solver.cc b/src/nufi_solver.cc
index deece5b..08c62a5 100644
--- a/src/nufi_solver.cc
+++ b/src/nufi_solver.cc
@@ -1,186 +1,187 @@
 #include "nufi/nufi_solver.h"
 
 #include <algorithm>
+#include <boost/qvm/mat_access.hpp>
 #include <cmath>
 #include <cstdlib>
 #include <deal.II/base/point.h>
 #include <deal.II/base/tensor.h>
 #include <deal.II/numerics/fe_field_function.h>
 
+#include <cstddef>
 #include <iostream>
 #include <memory>
 #include <ostream>
-#include <cstddef>
 #include <vector>
 
-#include "nufi/parameters.h"
-#include "nufi/save_results.h"
-#include "nufi/poisson_problem.h"
 #include "nufi/fields.h"
+#include "nufi/parameters.h"
+#include "nufi/poisson_problem.h"
+#include "nufi/save_results.h"
 #include "nufi/stopwatch.h"
 
 using namespace dealii;
 
-double NuFISolver::eval_ftilda(unsigned int n,
-                                      double x,
-                                      double u,
-                                      const PoissonProblem<1> &poisson) const
-{
-  if ( n == 0 ) return f0(x,u);
-  
+double
+NuFISolver::eval_ftilda(unsigned int n, double x, double u,
+                        const PoissonProblem<1> &poisson,
+                        const std::vector<Vector<double>> &phi_history) const {
+  if (n == 0)
+    return f0(x, u);
+
   double Ex;
 
   // We omit the initial half-step.
 
-  while ( --n )
-  {
-      x  = x - Parameters::DT *u;
-      Ex = -eval(x, poisson);
-      u  = u + Parameters::DT *Ex;
+  while (--n) {
+    x = x - Parameters::DT * u;
+    Ex = -eval(x, poisson, phi_history[n]);
+    u = u + Parameters::DT * Ex;
   }
 
   // The final half-step.
-  x -= Parameters::DT*u;
-  Ex = -eval(x, poisson);
-  u += 0.5*Parameters::DT*Ex;
+  x = x - Parameters::DT * u;
+  Ex = -eval(x, poisson, phi_history[n]);
+  u += 0.5 * Parameters::DT * Ex;
 
-  return f0(x,u);
+  return f0(x, u);
 }
 
-double NuFISolver::eval_f(unsigned int n,
-                                      double x,
-                                      double u,
-                                      const PoissonProblem<1> &poisson) const
-{
-  if ( n == 0 ) return f0(x,u);
-  
+double
+NuFISolver::eval_f(unsigned int n, double x, double u,
+                   const PoissonProblem<1> &poisson,
+                   const std::vector<Vector<double>> &phi_history) const {
+  if (n == 0)
+    return f0(x, u);
+
   double Ex;
 
   // Initial half-step.
-  Ex = -eval(x, poisson);
-  u += 0.5*Parameters::DT * Ex;
+  Ex = -eval(x, poisson, phi_history[n]);
+  u += 0.5 * Parameters::DT * Ex;
 
-  while ( --n )
-  {
-      x  = x - Parameters::DT *u;
-      Ex = -eval(x, poisson);
-      u  = u + Parameters::DT *Ex;
+  while (--n) {
+    x = x - Parameters::DT * u;
+    Ex = -eval(x, poisson, phi_history[n]);
+    u = u + Parameters::DT * Ex;
   }
 
   // The final half-step.
-  x -= Parameters::DT*u;
-  Ex = -eval(x, poisson);
-  u += 0.5*Parameters::DT*Ex;
+  x = x - Parameters::DT * u;
+  Ex = -eval(x, poisson, phi_history[n]);
+  u += 0.5 * Parameters::DT * Ex;
 
-  return f0(x,u);
+  return f0(x, u);
 }
 
-double NuFISolver::eval_rho(const unsigned int n,
-                          const double x,
-                          const PoissonProblem<1> &poisson,
-                          const unsigned int Nv) const
-{
-  const double dv = (Parameters::V_DOMAIN_RIGHT - Parameters::V_DOMAIN_LEFT) / Nv;
+double NuFISolver::eval_rho(const unsigned int n, const double x,
+                            const PoissonProblem<1> &poisson,
+                            const std::vector<Vector<double>> &phi_history,
+                            const unsigned int Nv) const {
+  const double dv =
+      (Parameters::V_DOMAIN_RIGHT - Parameters::V_DOMAIN_LEFT) / Nv;
   const double v_min = Parameters::V_DOMAIN_LEFT;
 
   double integral = 0.0;
 
-#pragma omp parallel for reduction (+ : integral)
-  for (unsigned int i = 0; i < Nv; ++i)
-    integral += eval_ftilda(n, x, v_min + i * dv, poisson);
-  
-  return 1.0 - integral*dv;
+  if (n == 0) {
+    for (unsigned int i = 0; i < Nv; ++i)
+      integral += f0(x, v_min + i * dv);
+  } else {
+#pragma omp parallel for reduction(+ : integral)
+    for (unsigned int i = 0; i < Nv; ++i)
+      integral += eval_ftilda(n, x, v_min + i * dv, poisson, phi_history);
+  }
+  return 1.0 - integral * dv;
 }
 
-void NuFISolver::run()
-{
+void NuFISolver::run() {
   std::cout << "Building E_sline\n\n";
 
   using std::abs;
   using std::max;
 
-  std::unique_ptr<double,decltype(std::free)*> rho { reinterpret_cast<double*>(std::aligned_alloc(64,sizeof(double)*Nx)), std::free };
+  std::unique_ptr<double, decltype(std::free) *> rho{
+      reinterpret_cast<double *>(std::aligned_alloc(64, sizeof(double) * Nx)),
+      std::free};
 
   std::vector<double> int_E_squared;
   int_E_squared.reserve(Nt);
 
-  if ( rho == nullptr ) throw std::bad_alloc {};
+  std::vector<Vector<double>> phi_history;
 
-  Gradient grad(x_min, x_max, Nx);
+  if (rho == nullptr)
+    throw std::bad_alloc{};
 
   double total_time = 0;
+  std::ofstream time_file("results/simulation_time.txt");
 
-  for (unsigned int it = 0; it < Nt; ++it)
-  {
-      stopwatch<double> timer;
-      
-      double time_elapsed_before = timer.elapsed();
+  for (unsigned int it = 0; it < Nt; ++it) {
+    stopwatch<double> timer;
 
-      std::cout << "Timestep " << it << " / " << Nt << " (simulation time = "<< it*Parameters::DT << ")"<< std::endl;
+    double time_elapsed_before = timer.elapsed();
 
-      // compute rho
+    std::cout << "Timestep " << it << " / " << Nt
+              << " (simulation time = " << it * Parameters::DT << ")"
+              << std::endl;
 
-      double dx = Parameters::CALC_DX;
-      
-      #pragma omp parallel for
-    	for(size_t i = 0; i<Nx; i++)
-    	{
-        double x = Parameters::X_DOMAIN_LEFT + i*dx;
-        double ith_rho = eval_rho(it, x, poisson, Parameters::NV); 
+    // compute rho
 
-        AssertThrow(std::isfinite(ith_rho), ExcMessage("NaN detected in rho"));
-    		rho.get()[i] =  ith_rho;
-    	}
+    double dx = Parameters::CALC_DX;
 
-      poisson.set_rhs_function(std::make_unique<ChargeDensity_NuFI<1>>(rho.get(), Nx));
-      poisson.solve_step();
+#pragma omp parallel for
+    for (size_t i = 0; i < Nx; i++) {
+      double x = Parameters::X_DOMAIN_LEFT + i * dx;
+      double ith_rho = eval_rho(it, x, poisson, phi_history, Parameters::NV);
 
-      std::vector<double> sampled_potential = poisson.sample_electric_potential(x_min, x_max, Nx); // Solution of FE
+      AssertThrow(std::isfinite(ith_rho), ExcMessage("NaN detected in rho"));
+      rho.get()[i] = ith_rho;
+    }
 
-      // These have been tested to be equivalent
-      // ////////////////////////////////////////////////
-      // std::vector<double> E_vals = grad.compute(sampled_potential); // vector grad of FE solution
-      // save_space_vector(E_vals, "electric", it);
-      // std::vector<double> E_vals_deal = poisson.sample_electric_field(x_min, x_max, Nx); // FE grad of soution
-      // save_space_vector(E_vals_deal, "electricdeal", it);
-      // ////////////////////////////////////////////////                                             
+    poisson.set_rhs_function(
+        std::make_unique<ChargeDensity_NuFI<1>>(rho.get(), Nx));
+    poisson.solve_step();
 
+    phi_history.push_back(poisson.get_solution());
 
-      // interpolate and save current field
+    std::vector<double> sampled_potential =
+        poisson.sample_electric_potential(x_min, x_max, Nx); // Solution of FE
 
-      std::vector<double> E_x(Nx,0.0) ;
-      #pragma omp parallel for
-      for(size_t ix=0; ix<Nx; ++ix)
-      {
-        E_x[ix] = -eval(Parameters::X_DOMAIN_LEFT+ix*dx, poisson);
+    double timer_elapsed = timer.elapsed();
+    double step_time = timer_elapsed - time_elapsed_before;
+    total_time += timer_elapsed;
+
+    std::cout << "step made in " << step_time << " seconds\n\n";
+    if (it % Parameters::PLOT_FREQUENCY == 0) {
+      std::cout << "Saving results...   ";
+      save_f(*this, it, poisson, phi_history, Parameters::PLOT_NX,
+             Parameters::NV, "results/ftilda_" + std::to_string(it) + ".dat");
+      save_rho(*this, it, poisson, phi_history, Parameters::PLOT_NX,
+               "results/rho_" + std::to_string(it) + ".dat");
+      // save_Efield(it, coeffs.get(), 128, "results/field_" +
+      // std::to_string(it) + ".dat");
+
+      std::vector<double> E_x(Nx, 0.0);
+#pragma omp parallel for
+      for (size_t ix = 0; ix < Nx; ++ix) {
+        E_x[ix] = -eval(Parameters::X_DOMAIN_LEFT + ix * dx, poisson,
+                        phi_history[it]);
       }
 
-      double timer_elapsed = timer.elapsed();
-      total_time += timer_elapsed;
-      
-      std::cout << "step made in "<< timer_elapsed-time_elapsed_before <<" seconds\n\n";
-      if (it % Parameters::PLOT_FREQUENCY == 0)
-      {
-        std::cout << "Saving results...   ";
-        save_f(*this, it, poisson, Parameters::PLOT_NX, Parameters::NV, "results/ftilda_" + std::to_string(it) + ".dat");
-        save_rho(*this, it, poisson, Parameters::PLOT_NX, "results/rho_" + std::to_string(it) + ".dat");
-        // save_Efield(it, coeffs.get(), 128, "results/field_" + std::to_string(it) + ".dat");
-        save_space_vector(E_x, "field", it);
+      save_space_vector(E_x, "field", it);
 
-        double int_val = 0.5 * integral_space_vector_squared(poisson);
-        int_E_squared.push_back(int_val);  
-        save_space_vector(int_E_squared, "electricint", it);
-        std::cout << "Time since start = "<< total_time<<"\n\n";
-      }
+      double int_val =
+          0.5 * integral_space_vector_squared(poisson, phi_history[it]);
+      int_E_squared.push_back(int_val);
+      save_space_vector(int_E_squared, "electricint", it);
+      std::cout << "Time since start = " << total_time << "\n\n";
+    }
   }
 
-  std::cout << "NuFI simulation finished in "<< total_time <<" seconds.\n";
+  std::cout << "NuFI simulation finished in " << total_time << " seconds.\n";
 }
 
-NuFISolver::NuFISolver()
-  : order(Parameters::FE_DEGREE),
-    poisson(order)
-{
+NuFISolver::NuFISolver() : order(Parameters::FE_DEGREE), poisson(order) {
   std::cout << "Initializing dealii Poisson Solver\n";
   poisson.initialize();
 }
diff --git a/src/save_results.cc b/src/save_results.cc
index 6d1b31d..3f32acd 100644
--- a/src/save_results.cc
+++ b/src/save_results.cc
@@ -1,22 +1,18 @@
 #include "nufi/save_results.h"
 
+#include "nufi/fields.h"
+#include "nufi/nufi_solver.h"
 #include "nufi/parameters.h"
+#include "nufi/poisson_problem.h"
 #include <fstream>
 #include <stdexcept>
 #include <string>
 #include <vector>
-#include "nufi/nufi_solver.h"
-#include "nufi/poisson_problem.h"
-#include "nufi/fields.h"
 
-
-void save_f( const NuFISolver &solver,
-    unsigned int n,
-                                    const PoissonProblem<1> &poisson,
-                                    unsigned int Nx_out,
-                                    unsigned int Nv_out,
-                                    const std::string &filename)
-{
+void save_f(const NuFISolver &solver, unsigned int n,
+            const PoissonProblem<1> &poisson,
+            const std::vector<Vector<double>> &phi_history, unsigned int Nx_out,
+            unsigned int Nv_out, const std::string &filename) {
   std::ofstream file(filename);
 
   double xmin = Parameters::X_DOMAIN_LEFT;
@@ -32,34 +28,30 @@ void save_f( const NuFISolver &solver,
   file << xmin << " " << xmax << "\n";
   file << vmin << " " << vmax << "\n";
 
-  for (unsigned int i = 0; i < Nx_out; ++i)
-  {
-      double x = xmin + (i + 0.5)*dx;
+  for (unsigned int i = 0; i < Nx_out; ++i) {
+    double x = xmin + (i + 0.5) * dx;
 
-      for (unsigned int j = 0; j < Nv_out; ++j)
-      {
-          double v = vmin + (j + 0.5)*dv;
+    for (unsigned int j = 0; j < Nv_out; ++j) {
+      double v = vmin + (j + 0.5) * dv;
 
-          double val = solver.eval_f(n, x, v, poisson);
+      double val = solver.eval_f(n, x, v, poisson, phi_history);
 
-          file << val;
+      file << val;
 
-          if (j < Nv_out - 1)
-              file << " ";
-      }
+      if (j < Nv_out - 1)
+        file << " ";
+    }
 
-      file << "\n";
+    file << "\n";
   }
 
   file.close();
 }
 
-void save_rho(const NuFISolver &solver,
-    unsigned int n,
-                                    const PoissonProblem<1> &poisson,
-                                    unsigned int Nx_out,
-                                    const std::string &filename)
-{
+void save_rho(const NuFISolver &solver, unsigned int n,
+              const PoissonProblem<1> &poisson,
+              const std::vector<Vector<double>> &phi_history,
+              unsigned int Nx_out, const std::string &filename) {
   std::ofstream file(filename);
 
   double xmin = Parameters::X_DOMAIN_LEFT;
@@ -69,20 +61,18 @@ void save_rho(const NuFISolver &solver,
   file << Nx_out << "\n";
   file << xmin << " " << xmax << "\n";
 
-  for (unsigned int i = 0; i < Nx_out; ++i, xmin += dx)
-  {
-      double val = solver.eval_rho(n, xmin, poisson);
-      file << val;
-      file << "\n";
+  for (unsigned int i = 0; i < Nx_out; ++i, xmin += dx) {
+    double val = solver.eval_rho(n, xmin, poisson, phi_history);
+    file << val;
+    file << "\n";
   }
   file.close();
 }
 
-void save_Efield([[maybe_unused]]unsigned int n,
-                                    const PoissonProblem<1> &poisson,
-                                    unsigned int Nx_out,
-                                    const std::string &filename)
-{
+void save_Efield([[maybe_unused]] unsigned int n,
+                 const PoissonProblem<1> &poisson,
+                 const std::vector<Vector<double>> &phi_history,
+                 unsigned int Nx_out, const std::string &filename) {
   std::ofstream file(filename);
 
   double xmin = Parameters::X_DOMAIN_LEFT;
@@ -94,23 +84,25 @@ void save_Efield([[maybe_unused]]unsigned int n,
   file << Nx_out << "\n";
   file << xmin << " " << xmax << "\n";
 
-  for (unsigned int i = 0; i < Nx_out; ++i, xmin += dx)
-  {
-      double val = -eval(xmin, poisson);
-      file << val;
-      file << "\n";
+  for (unsigned int i = 0; i < Nx_out; ++i, xmin += dx) {
+    double val = -eval(xmin, poisson, phi_history[n]);
+    file << val;
+    file << "\n";
   }
   file.close();
 }
 
-void save_space_vector(const std::vector<double>& vals, const std::string& filename, size_t it)
-{
-    std::ofstream file("results/" + filename + "_" + std::to_string(it) + ".dat");
+void save_space_vector(const std::vector<double> &vals,
+                       const std::string &filename, size_t it) {
+  std::ofstream file("results/" + filename + "_" + std::to_string(it) + ".dat");
 
-    if (!file) throw std::runtime_error("failed to start file in results/");
+  if (!file)
+    throw std::runtime_error("failed to start file in results/");
 
-    file << vals.size() << "\n";
-    file << Parameters::X_DOMAIN_LEFT << " " << Parameters::X_DOMAIN_RIGHT << "\n";
-    file << std::fixed << std::setprecision(8);
-    for (double val : vals) file << val << "\n";
+  file << vals.size() << "\n";
+  file << Parameters::X_DOMAIN_LEFT << " " << Parameters::X_DOMAIN_RIGHT
+       << "\n";
+  file << std::fixed << std::setprecision(8);
+  for (double val : vals)
+    file << val << "\n";
 }