diff --git a/libnufi_lib.a b/libnufi_lib.a
index dd03ba6..a9045c0 100644
Binary files a/libnufi_lib.a and b/libnufi_lib.a differ
diff --git a/nufi/cells.h b/nufi/cells.h
index 36c61a2..95a8ac2 100644
--- a/nufi/cells.h
+++ b/nufi/cells.h
@@ -2,6 +2,7 @@
 #define CELLS_H
 
 #include <algorithm>
+#include <boost/geometry/geometries/concepts/point_concept.hpp>
 #include <deal.II/base/geometry_info.h>
 #include <deal.II/base/point.h>
 #include <deal.II/dofs/dof_handler.h>
@@ -32,8 +33,11 @@ public:
                const Triangulation<dim> &triangulation);
   CellLocation<dim> locate(const Point<dim> &p) const;
 
+  const std::vector<Point<dim>> &get_cell_centers() const;
+
 private:
   std::vector<CellInfo<dim>> cells;
+  std::vector<Point<dim>> cell_centers;
 };
 
 template <int dim>
@@ -59,6 +63,17 @@ void CellLocator<dim>::rebuild(const DoFHandler<dim> &dof_handler,
             [](const CellInfo<dim> &a, const CellInfo<dim> &b) {
               return a.lower[0] < b.lower[0];
             });
+
+  cell_centers.clear();
+  cell_centers.reserve(cells.size());
+
+  for (const auto &cell : cells) {
+    Point<dim> center;
+    for (unsigned int d = 0; d < dim; ++d)
+      center[d] = 0.5 * (cell.lower[d] + cell.upper[d]);
+
+    cell_centers.push_back(center);
+  }
 }
 
 template <int dim>
@@ -93,4 +108,9 @@ CellLocation<dim> CellLocator<dim>::locate(const Point<dim> &p) const {
   return location;
 }
 
+template <int dim>
+const std::vector<Point<dim>> &CellLocator<dim>::get_cell_centers() const {
+  return cell_centers;
+}
+
 #endif // !CELLS_H
diff --git a/nufi/fields.h b/nufi/fields.h
index c8951be..46e58f6 100644
--- a/nufi/fields.h
+++ b/nufi/fields.h
@@ -61,15 +61,17 @@ using namespace dealii;
 // }
 inline std::vector<double> make_x_eval(size_t Nx) {
   std::vector<double> x_eval(Nx);
+  const double dx = Parameters::LX / Nx;
   for (size_t i = 0; i < Nx; ++i)
-    x_eval[i] = Parameters::X_DOMAIN_LEFT + i * Parameters::CALC_DX;
+    x_eval[i] = Parameters::X_DOMAIN_LEFT + i * dx;
   return x_eval;
 }
 
 inline void reset_x_eval(std::vector<double> &x_vals) {
   const size_t Nx = x_vals.size();
+  const double dx = Parameters::LX / Nx;
   for (size_t i = 0; i < Nx; ++i)
-    x_vals[i] = Parameters::X_DOMAIN_LEFT + i * Parameters::CALC_DX;
+    x_vals[i] = Parameters::X_DOMAIN_LEFT + i * dx;
 };
 
 inline double f0(const double x, const double v,
@@ -132,4 +134,15 @@ inline double integral_space_vector_squared(const PoissonProblem<1> &poisson,
   return integral * dx;
 };
 
+inline std::vector<double>
+Point_vector_to_double_vector(const std::vector<Point<1>> &Points) {
+  const size_t n_points = Points.size();
+  std::vector<double> vector(n_points);
+
+  for (size_t i = 0; i < n_points; ++i)
+    vector[i] = Points[i][0];
+
+  return vector;
+}
+
 #endif
diff --git a/nufi/poisson_problem.h b/nufi/poisson_problem.h
index 2db18db..641cdb2 100644
--- a/nufi/poisson_problem.h
+++ b/nufi/poisson_problem.h
@@ -70,7 +70,10 @@ public:
                               std::vector<Vector<double>> &solution_history);
   void run();
 
+  unsigned int get_rhs_size();
+  unsigned int get_dof_size();
   void set_rhs_function(std::function<double(const Point<dim> &)> f);
+  void set_rhs(const Vector<double> &new_rhs) { rhs = new_rhs; }
 
   const Vector<double> &get_solution() const { return solution; }
   const MappingQ<dim> &get_mapping() const { return mapping; }
@@ -89,6 +92,7 @@ public:
 
   Triangulation<dim> triangulation;
   DoFHandler<dim> dof_handler;
+  CellLocator<dim> cell_locator;
 
 private:
   void create_mesh();
@@ -109,12 +113,10 @@ private:
   Vector<double> system_rhs;
 
   std::function<double(const Point<dim> &)> rhs_function;
+  Vector<double> rhs;
 
   MappingQ<dim> mapping;
 
-  CellLocator<dim> cell_locator;
-  // std::vector<typename DoFHandler<dim>::active_cell_iterator> active_cells;
-
   mutable std::vector<double> local_solution_buffer;
   mutable std::unique_ptr<FEPointEvaluation<dim, dim>> evaluator;
 };
@@ -123,12 +125,22 @@ private:
 // Utilities
 //====//====//
 
-template <int dim>
-void PoissonProblem<dim>::set_rhs_function(
-    std::function<double(const Point<dim> &)> f) {
-  rhs_function = std::move(f);
+template <int dim> unsigned int PoissonProblem<dim>::get_rhs_size() {
+
+  QGauss<dim> quadrature_formula(fe.degree + 1);
+  return quadrature_formula.size();
 }
 
+template <int dim> unsigned int PoissonProblem<dim>::get_dof_size() {
+  return dof_handler.n_dofs();
+}
+
+// template <int dim>
+// void PoissonProblem<dim>::set_rhs_function(
+//     std::function<double(const Point<dim> &)> f) {
+//   rhs_function = std::move(f);
+// }
+
 template <int dim>
 PoissonProblem<dim>::PoissonProblem(unsigned int degree)
     : triangulation(Triangulation<dim>::limit_level_difference_at_vertices),
@@ -358,15 +370,19 @@ 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);
 
+  std::vector<double> rhs_values(quadrature_formula.size());
+
   for (const auto &cell : dof_handler.active_cell_iterators()) {
     fe_values.reinit(cell);
 
     cell_matrix = 0;
     cell_rhs = 0;
 
+    fe_values.get_function_values(rhs, rhs_values);
+
     for (const auto q : fe_values.quadrature_point_indices()) {
-      const double rho = rhs_function(
-          fe_values.quadrature_point(q)); // Eval rhs_function at q points
+      // 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())
@@ -374,7 +390,9 @@ template <int dim> void PoissonProblem<dim>::assemble_system() {
                                fe_values.shape_grad(j, q) * fe_values.JxW(q);
 
       for (const unsigned int i : fe_values.dof_indices())
-        cell_rhs(i) += fe_values.shape_value(i, q) * rho * fe_values.JxW(q);
+        // cell_rhs(i) += fe_values.shape_value(i, q) * rho * fe_values.JxW(q);
+        cell_rhs(i) +=
+            fe_values.shape_value(i, q) * rhs_values[q] * fe_values.JxW(q);
     }
 
     cell->get_dof_indices(local_dof_indices);
@@ -433,7 +451,8 @@ void PoissonProblem<dim>::coarse_and_refine_grid(
 template <int dim> void PoissonProblem<dim>::solve() {
 
   SolverControl solver_control(Parameters::CONVERGENCE_ITERATIONS,
-                               Parameters::CONVERGENCE_LIMIT);
+                               Parameters::CONVERGENCE_LIMIT *
+                                   system_rhs.l2_norm());
   SolverCG<Vector<double>> solver(solver_control);
 
   // PreconditionSSOR<SparseMatrix<double>> preconditioner;
diff --git a/src/nufi_solver.cc b/src/nufi_solver.cc
index be078a4..0d6131f 100644
--- a/src/nufi_solver.cc
+++ b/src/nufi_solver.cc
@@ -186,8 +186,8 @@ void NuFISolver::run() {
             << "plot_time"
             << "\n";
 
-  const double x_min = Parameters::X_DOMAIN_LEFT;
-  double dx = Parameters::CALC_DX;
+  [[maybe_unused]] const double x_min = Parameters::X_DOMAIN_LEFT;
+  [[maybe_unused]] double dx = Parameters::CALC_DX;
 
   //====//====//
   // Time loop//
@@ -222,21 +222,15 @@ void NuFISolver::run() {
 
     double compute_start = timer.elapsed();
     // compute rho
-    std::vector<double> x_eval = make_x_eval(Nx);
-    std::vector<double> tmp_rho =
+    std::vector<double> x_eval = make_x_eval(poisson.get_dof_size());
+    std::vector<double> rho_values =
         eval_rho(it, x_eval, poisson, phi_history, Parameters::NV);
 
-    for (size_t i = 0; i < Nx; i++) {
-      AssertThrow(std::isfinite(tmp_rho[i]), ExcMessage("NaN detected in rho"));
-      rho.get()[i] = tmp_rho[i];
-    }
+    Vector<double> rhs(x_eval.size());
+    for (unsigned int i = 0; i < rhs.size(); ++i)
+      rhs[i] = rho_values[i];
 
-    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.set_rhs(rhs);
 
     poisson.solve_step();
     phi_history.push_back(poisson.get_solution());
@@ -248,7 +242,8 @@ void NuFISolver::run() {
       poisson.coarse_and_refine_grid(it, phi_history);
       refine_time = timer.elapsed() - refine_start;
       std::cout << "Refinement step done in "
-                << std::to_string(std::floor(refine_time)) << "[s]" << "\n";
+                << std::to_string(std::round(std::floor(refine_time))) << "[s]"
+                << "\n";
     }
 
     double timer_elapsed = timer.elapsed();
@@ -266,6 +261,7 @@ void NuFISolver::run() {
       // std::to_string(it) + ".dat");
 
       std::vector<double> x_eval_Ex = make_x_eval(Parameters::PLOT_NX);
+      std::vector<double> tmp_rho(x_eval_Ex.size());
       tmp_rho = eval(x_eval_Ex, poisson, phi_history[it]);
 
       std::vector<double> E_x(Parameters::PLOT_NX);