diff --git a/libnufi_lib.a b/libnufi_lib.a
index 0b9082d..fb38e5b 100644
Binary files a/libnufi_lib.a and b/libnufi_lib.a differ
diff --git a/nufi/fields.h b/nufi/fields.h
index 84ed5a7..97aa7cf 100644
--- a/nufi/fields.h
+++ b/nufi/fields.h
@@ -2,7 +2,7 @@
 #define FIELDS_H
 
 #include "nufi/parameters.h"
-#include "poisson_problem.h"
+#include "nufi/poisson_problem.h"
 #include <cmath>
 #include <deal.II/base/function.h>
 #include <deal.II/base/point.h>
@@ -68,27 +68,14 @@ inline double f0(const double x, const double v,
   return prefactor * gaussian;
 }
 
-inline double compute_rho(const double x,
-                          const unsigned int Nv = Parameters::NV) {
-  const double dv =
-      (Parameters::V_DOMAIN_RIGHT - Parameters::V_DOMAIN_LEFT) / Nv;
-
-  double integral = 0.0;
-
-  for (unsigned int i = 0; i < Nv; ++i) {
-    const double v = Parameters::V_DOMAIN_LEFT + (i + 0.5) * dv;
-    integral += f0(x, v) * dv;
-  }
-
-  return 1.0 - integral;
-}
-
+// wrapper for eval_point() { VectorTools::point_values() }
 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);
+  x = x - Parameters::LX * std::floor(x * Parameters::LX_INV); // in domain
+
   return eval_point<1>(poisson.get_mapping(), poisson.get_dof_handler(),
                        solution, Point<1>(x));
 }
@@ -122,58 +109,4 @@ inline double integral_space_vector_squared(const PoissonProblem<1> &poisson,
   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_;
-// };
-
-template <int dim>
-class ChargeDensity : public Function<dim> // only uses f0
-{
-public:
-  ChargeDensity(double eps, double k, unsigned int Nv)
-      : Function<dim>(1), eps(eps), k(k), Nv(Nv) {}
-
-  virtual double
-  value(const Point<dim> &p,
-        [[maybe_unused]] const unsigned int component = 0) const override {
-    return compute_rho(p[0], Nv);
-  }
-
-private:
-  const double eps;
-  const double k;
-  const unsigned int Nv;
-};
-
 #endif
diff --git a/nufi/nufi_solver.h b/nufi/nufi_solver.h
index 9af3813..95fc301 100644
--- a/nufi/nufi_solver.h
+++ b/nufi/nufi_solver.h
@@ -5,7 +5,7 @@
 #include <cmath>
 #include <deal.II/base/point.h>
 #include <deal.II/base/tensor.h>
-#include <deal.II/numerics/fe_field_function.h>
+#include <deal.II/numerics/vector_tools.h>
 #include <vector>
 
 #include "nufi/fields.h" //dont remove
@@ -19,9 +19,10 @@ public:
   NuFISolver();
 
   void run();
-  double eval_rho(unsigned int n, double x, const PoissonProblem<1> &poisson,
+  double eval_rho(unsigned int n, const double x,
+                  const PoissonProblem<1> &poisson,
                   const std::vector<Vector<double>> &phi_history,
-                  unsigned int Nv = Parameters::NV) const;
+                  const 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;
@@ -44,32 +45,4 @@ private:
 
   PoissonProblem<1> poisson;
 };
-
-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];
-
-    // 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));
-
-    // periodic wrap
-    i = (i % Nx + Nx) % Nx;
-
-    return rho[i];
-  }
-
-private:
-  const double *rho;
-  const unsigned int Nx;
-};
-
 #endif
diff --git a/nufi/parameters.h b/nufi/parameters.h
index a457124..3b533f8 100644
--- a/nufi/parameters.h
+++ b/nufi/parameters.h
@@ -32,7 +32,7 @@ constexpr double WAVE_NR = 0.5;
 constexpr double F0_FACTOR = 0.39894228040143267793994; // 1/sqrt(2pi)
 
 // NUFI options
-constexpr double DT = 1. / 10.;
+constexpr double DT = 1. / 8.;
 constexpr unsigned int TMAX = 100;
 
 // Plotting options
diff --git a/nufi/poisson_problem.h b/nufi/poisson_problem.h
index cae16a3..573ebf2 100644
--- a/nufi/poisson_problem.h
+++ b/nufi/poisson_problem.h
@@ -33,12 +33,12 @@
 #include <deal.II/fe/fe_values.h>
 
 #include <deal.II/numerics/data_out.h>
-#include <deal.II/numerics/fe_field_function.h>
 #include <deal.II/numerics/matrix_tools.h>
 #include <deal.II/numerics/vector_tools.h>
 
 #include <deal.II/numerics/vector_tools_evaluate.h>
 #include <deal.II/numerics/vector_tools_interpolate.h>
+#include <functional>
 #include <memory>
 #include <string>
 #include <utility>
@@ -58,11 +58,11 @@ public:
   void solve_step();
   void run();
 
-  void set_rhs_function(std::unique_ptr<Function<dim>> rhs_function);
+  void set_rhs_function(std::function<double(const Point<dim> &)> f);
 
   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; }
+  const DoFHandler<dim> &get_dof_handler() const { return dof_handler; }
 
   std::vector<double> sample_electric_field(double x_min, double x_max,
                                             unsigned int Nx);
@@ -87,18 +87,17 @@ private:
   Vector<double> solution; // phi
   Vector<double> system_rhs;
 
-  std::unique_ptr<const Function<dim>> rhs_function;
+  std::function<double(const Point<dim> &)> rhs_function;
 
   MappingQ<dim> mapping;
-
-  // std::unique_ptr<Functions::FEFieldFunction<dim>> fe_field_function;
 };
 
 // Utilities
 
 template <int dim>
-void PoissonProblem<dim>::set_rhs_function(std::unique_ptr<Function<dim>> rhs) {
-  rhs_function = std::move(rhs);
+void PoissonProblem<dim>::set_rhs_function(
+    std::function<double(const Point<dim> &)> f) {
+  rhs_function = std::move(f);
 }
 
 template <int dim>
@@ -170,71 +169,6 @@ PoissonProblem<dim>::sample_electric_potential(double x_min, double x_max,
   return values;
 }
 
-// // by GPT to re-re-re-check
-// 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 DoFHandler<dim> &dof_handler,
-//         const MappingQ<dim> &mapping)
-// {
-//     AssertIndexRange(n, solutions.size());
-//     Assert(points.size() == cell_indices.size(),
-//            ExcMessage("points and cell_indices must have same size"));
-//
-//     const Vector<double> &solution = solutions[n];
-//
-//     std::vector<double> result(points.size());
-//
-//     // Group points by cell (required for FEPointEvaluation efficiency)
-//     std::map<unsigned int, std::vector<unsigned int>> cell_to_point_ids;
-//
-//     for (unsigned int i = 0; i < points.size(); ++i)
-//         cell_to_point_ids[cell_indices[i]].push_back(i);
-//
-//     FEPointEvaluation<1, dim> evaluator(mapping,
-//                                         dof_handler.get_fe(),
-//                                         update_values);
-//
-//     std::vector<Point<dim>> cell_points;
-//     Vector<double> local_dofs(dof_handler.get_fe().dofs_per_cell);
-//
-//     for (const auto &entry : cell_to_point_ids)
-//     {
-//         const unsigned int cell_id = entry.first;
-//         const auto &point_ids = entry.second;
-//
-//         // these two lines bellow assume some order not sure how or why
-//         auto cell = dof_handler.begin_active();
-//         std::advance(cell, cell_id);
-//
-//         // extract points belonging to this cell
-//         cell_points.clear();
-//         cell_points.reserve(point_ids.size());
-//
-//         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());
-//         cell->get_dof_indices(indices);
-//
-//         for (unsigned int i=0;i<indices.size();++i)
-//           local_dofs[i] = solution[indices[i]];
-//
-//         // initialize evaluator on this cell
-//         evaluator.reinit(cell, cell_points);
-//
-//         evaluator.evaluate(local_dofs, EvaluationFlags::values);
-//
-//         for (unsigned int k = 0; k < point_ids.size(); ++k)
-//             result[point_ids[k]] = evaluator.get_value(k);
-//     }
-//
-//     return result;
-// }
-
 template <int dim>
 double eval_point(const Mapping<dim> &mapping,
                   const DoFHandler<dim> &dof_handler,
@@ -299,90 +233,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);
 }
-/* (Mine)
-template <int dim>
-void PoissonProblem<dim>::assemble_system()
-{
-
-  system_matrix = 0;
-  system_rhs = 0;
-
-  QGauss<dim>  quadrature_formula(fe.degree + 1);
-  FEValues<dim> fe_values(fe, quadrature_formula,
-                          update_values |
-                          update_gradients |
-                          update_quadrature_points |
-                          update_JxW_values);
-
-  const unsigned int dofs_per_cell = fe.n_dofs_per_cell();
-
-  FullMatrix<double> cell_matrix(dofs_per_cell, dofs_per_cell);
-  Vector<double>     cell_rhs(dofs_per_cell);
-  std::vector<types::global_dof_index> local_dof_indices(dofs_per_cell);
-
-  Assert(rhs_function != nullptr, ExcMessage("RHS function not set"));
-
-  for (const auto &cell : dof_handler.active_cell_iterators())
-  {
-    fe_values.reinit(cell);
-
-    cell_matrix = 0;
-    cell_rhs    = 0;
-
-    for (const auto q : fe_values.quadrature_point_indices())
-    {
-      const double rho = rhs_function->value(fe_values.quadrature_point(q));
-
-      for (const unsigned int i : fe_values.dof_indices())
-        for (const unsigned int j : fe_values.dof_indices())
-          cell_matrix(i, j) +=
-            (fe_values.shape_grad(i, q) * // grad phi_i(x_q)
-             fe_values.shape_grad(j, q) * // grad phi_j(x_q)
-             fe_values.JxW(q));           // dx
-
-      for (const unsigned int i : fe_values.dof_indices())
-        cell_rhs(i) += (fe_values.shape_value(i, q) * // phi_i(x_q)
-                        rho *                         // f(x_q)
-                        fe_values.JxW(q));            // dx
-
-
-    }
-
-    cell->get_dof_indices(local_dof_indices);
-    constraints.distribute_local_to_global(cell_matrix,
-                                           cell_rhs,
-                                           local_dof_indices,
-                                           system_matrix,
-                                           system_rhs);
-    for (const unsigned int i : fe_values.dof_indices())
-      for (const unsigned int j : fe_values.dof_indices())
-        system_matrix.add(local_dof_indices[i],
-                          local_dof_indices[j],
-                          cell_matrix(i, j));
-
-    for (const unsigned int i : fe_values.dof_indices())
-      system_rhs(local_dof_indices[i]) += cell_rhs(i);
-
-  }
-  std::map<types::global_dof_index, double> boundary_values;
-  // VectorTools::interpolate_boundary_values(dof_handler,
-  //                                          types::boundary_id(0),
-  //                                          Functions::ZeroFunction<1>(),
-  //                                          boundary_values);
-  MatrixTools::apply_boundary_values(boundary_values,
-                                     system_matrix,
-                                     solution,
-                                     system_rhs);
-}
- */
-
-// Paul's, mine's above
 
+// Paul
 template <int dim> void PoissonProblem<dim>::assemble_system() {
   Assert(system_matrix.m() == dof_handler.n_dofs(),
          ExcMessage("Matrix not initialized correctly"));
@@ -400,8 +253,6 @@ template <int dim> void PoissonProblem<dim>::assemble_system() {
   Vector<double> cell_rhs(dofs_per_cell);
   std::vector<types::global_dof_index> local_dof_indices(dofs_per_cell);
 
-  Assert(rhs_function != nullptr, ExcMessage("RHS function not set"));
-
   for (const auto &cell : dof_handler.active_cell_iterators()) {
     fe_values.reinit(cell);
 
@@ -409,7 +260,8 @@ template <int dim> void PoissonProblem<dim>::assemble_system() {
     cell_rhs = 0;
 
     for (const auto q : fe_values.quadrature_point_indices()) {
-      const double rho = rhs_function->value(fe_values.quadrature_point(q));
+      const double rho = rhs_function(
+          fe_values.quadrature_point(q)); // Eval rhs_function at q points
 
       for (const unsigned int i : fe_values.dof_indices())
         for (const unsigned int j : fe_values.dof_indices())
@@ -437,6 +289,7 @@ template <int dim> void PoissonProblem<dim>::solve() {
   // preconditioner.initialize(system_matrix, 1.2);
 
   // solver.solve(system_matrix, solution, system_rhs, preconditioner);
+
   solver.solve(system_matrix, solution, system_rhs, PreconditionIdentity());
   constraints.distribute(solution);
 }
diff --git a/rho_E.mp4 b/rho_E.mp4
deleted file mode 100644
index 2730435..0000000
Binary files a/rho_E.mp4 and /dev/null differ
diff --git a/src/nufi_solver.cc b/src/nufi_solver.cc
index 08c62a5..ac90569 100644
--- a/src/nufi_solver.cc
+++ b/src/nufi_solver.cc
@@ -6,7 +6,7 @@
 #include <cstdlib>
 #include <deal.II/base/point.h>
 #include <deal.II/base/tensor.h>
-#include <deal.II/numerics/fe_field_function.h>
+#include <deal.II/numerics/vector_tools.h>
 
 #include <cstddef>
 #include <iostream>
@@ -74,7 +74,7 @@ NuFISolver::eval_f(unsigned int n, double x, double u,
   return f0(x, u);
 }
 
-double NuFISolver::eval_rho(const unsigned int n, const double x,
+double NuFISolver::eval_rho(unsigned int n, const double x,
                             const PoissonProblem<1> &poisson,
                             const std::vector<Vector<double>> &phi_history,
                             const unsigned int Nv) const {
@@ -84,14 +84,9 @@ double NuFISolver::eval_rho(const unsigned int n, const double x,
 
   double integral = 0.0;
 
-  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);
-  }
+  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;
 }
 
@@ -115,6 +110,10 @@ void NuFISolver::run() {
 
   double total_time = 0;
   std::ofstream time_file("results/simulation_time.txt");
+  time_file << "# it step_time total_time" << "\n";
+
+  const double x_min = Parameters::X_DOMAIN_LEFT;
+  double dx = Parameters::CALC_DX;
 
   for (unsigned int it = 0; it < Nt; ++it) {
     stopwatch<double> timer;
@@ -127,8 +126,6 @@ void NuFISolver::run() {
 
     // compute rho
 
-    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;
@@ -138,8 +135,13 @@ void NuFISolver::run() {
       rho.get()[i] = ith_rho;
     }
 
-    poisson.set_rhs_function(
-        std::make_unique<ChargeDensity_NuFI<1>>(rho.get(), Nx));
+    poisson.set_rhs_function([&rho, x_min, dx, Nx = Nx](const Point<1> &p) {
+      double x = p[0];
+      int i = static_cast<int>(std::floor((x - x_min) / dx));
+      i = (i % Nx + Nx) % Nx;
+      return rho.get()[i];
+    });
+
     poisson.solve_step();
 
     phi_history.push_back(poisson.get_solution());
@@ -151,6 +153,8 @@ void NuFISolver::run() {
     double step_time = timer_elapsed - time_elapsed_before;
     total_time += timer_elapsed;
 
+    time_file << it << " " << step_time << " " << total_time << "\n";
+
     std::cout << "step made in " << step_time << " seconds\n\n";
     if (it % Parameters::PLOT_FREQUENCY == 0) {
       std::cout << "Saving results...   ";
@@ -164,8 +168,7 @@ void NuFISolver::run() {
       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]);
+        E_x[ix] = -eval(x_min + ix * dx, poisson, phi_history[it]);
       }
 
       save_space_vector(E_x, "field", it);