diff --git a/libnufi_lib.a b/libnufi_lib.a
index 5703b05..031c6a7 100644
Binary files a/libnufi_lib.a and b/libnufi_lib.a differ
diff --git a/nufi/nufi_solver.h b/nufi/nufi_solver.h
index f9a9d18..32de460 100644
--- a/nufi/nufi_solver.h
+++ b/nufi/nufi_solver.h
@@ -34,6 +34,12 @@ public:
          const std::vector<GridStructure<1>> &grid_struct,
          const std::vector<SolutionSnapshot<1>> &phi_history) const;
 
+  std::vector<double>
+  eval_rho_points(unsigned int n, const std::vector<Point<1>> &points,
+                  const std::vector<GridStructure<1>> &grid_struct,
+                  const std::vector<SolutionSnapshot<1>> &phi_history,
+                  const unsigned int Nv) const;
+
 private:
   unsigned int Nt = std::floor(Parameters::TMAX / Parameters::DT);
   unsigned int Nx = Parameters::CALC_NX;
diff --git a/nufi/poisson_problem.h b/nufi/poisson_problem.h
index 4a2517f..35a4404 100644
--- a/nufi/poisson_problem.h
+++ b/nufi/poisson_problem.h
@@ -73,8 +73,9 @@ public:
 
   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; }
+  void set_rhs_function(
+      std::function<std::vector<double>(const std::vector<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; }
@@ -116,8 +117,9 @@ private:
   Vector<double> solution; // phi
   Vector<double> system_rhs;
 
-  std::function<double(const Point<dim> &)> rhs_function;
-  Vector<double> rhs;
+  std::function<std::vector<double>(const std::vector<Point<dim>> &)>
+      rhs_function;
+  // Vector<double> rhs;
 
   MappingQ<dim> mapping;
 
@@ -139,11 +141,11 @@ 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>
+void PoissonProblem<dim>::set_rhs_function(
+    std::function<std::vector<double>(const std::vector<Point<dim>> &)> f) {
+  rhs_function = std::move(f);
+}
 
 template <int dim>
 PoissonProblem<dim>::PoissonProblem(unsigned int degree)
@@ -331,14 +333,14 @@ template <int dim> void PoissonProblem<dim>::setup_system() {
   // update_gradients);
 }
 
-// Paul
 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;
 
-  QGauss<dim> quadrature_formula(fe.degree + 1);
+  const QGauss<dim> quadrature_formula(fe.degree + 1);
   FEValues<dim> fe_values(fe, quadrature_formula,
                           update_values | update_gradients |
                               update_quadrature_points | update_JxW_values);
@@ -347,32 +349,49 @@ template <int dim> void PoissonProblem<dim>::assemble_system() {
 
   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);
 
-  std::vector<double> rhs_values(quadrature_formula.size());
+  // Eval rhs_function only once for all quadrature points
+  const unsigned int n_q_points = quadrature_formula.size();
+  std::vector<Point<dim>> all_q_points;
+  all_q_points.reserve(triangulation.n_active_cells() * n_q_points);
 
+  for (const auto &cell : dof_handler.active_cell_iterators()) {
+    fe_values.reinit(cell);
+    const auto &q_points = fe_values.get_quadrature_points();
+    all_q_points.insert(all_q_points.end(), q_points.begin(), q_points.end());
+  }
+
+  Assert(rhs_function,
+         ExcMessage("Poisson RHS function has not been initialized."));
+
+  std::vector<double> all_rho = rhs_function(all_q_points);
+
+  Assert(all_rho.size() == all_q_points.size(),
+         ExcMessage("rhs_function returned wrong size"));
+
+  // assemble system
+  unsigned int q_offset = 0;
   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 (size_t q = 0; q < n_q_points; ++q) {
 
-    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
-
-      for (const unsigned int i : fe_values.dof_indices())
-        for (const unsigned int j : fe_values.dof_indices())
+      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) *
                                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) * rhs_values[q] * fe_values.JxW(q);
+        cell_rhs(i) += fe_values.shape_value(i, q) * all_rho[q_offset + q] *
+                       fe_values.JxW(q);
+      }
     }
+    q_offset += n_q_points;
 
     cell->get_dof_indices(local_dof_indices);
 
diff --git a/src/nufi_solver.cc b/src/nufi_solver.cc
index 84f68fe..be03cfb 100644
--- a/src/nufi_solver.cc
+++ b/src/nufi_solver.cc
@@ -158,6 +158,15 @@ NuFISolver::eval_rho(unsigned int n, std::vector<double> &X,
   return integral;
 }
 
+std::vector<double>
+NuFISolver::eval_rho_points(unsigned int n, const std::vector<Point<1>> &points,
+                            const std::vector<GridStructure<1>> &grid_struct,
+                            const std::vector<SolutionSnapshot<1>> &phi_history,
+                            const unsigned int Nv) const {
+  std::vector<double> point_vector = Point_vector_to_double_vector(points);
+  return NuFISolver::eval_rho(n, point_vector, grid_struct, phi_history, Nv);
+}
+
 void NuFISolver::run() {
 
   //====//====//
@@ -167,12 +176,13 @@ void NuFISolver::run() {
   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};
-
-  if (rho == nullptr)
-    throw std::bad_alloc{};
+  // std::unique_ptr<double, decltype(std::free) *> rho{
+  //     reinterpret_cast<double *>(std::aligned_alloc(64, sizeof(double) *
+  //     Nx)),
+  // std::free};
+  //
+  // if (rho == nullptr)
+  //   throw std::bad_alloc{};
 
   std::vector<double> int_E_squared;
   int_E_squared.reserve(Nt);
@@ -234,15 +244,15 @@ void NuFISolver::run() {
 
     double compute_start = timer.elapsed();
     // compute rho
-    std::vector<double> x_eval = make_x_eval(poisson.get_dof_size());
-    std::vector<double> rho_values =
-        eval_rho(it, x_eval, grid_versions, phi_history, Parameters::NV);
+    //
+    poisson.set_rhs_function([&](const std::vector<Point<1>> &points) {
+      std::vector<double> x(points.size());
 
-    Vector<double> rhs(x_eval.size());
-    for (unsigned int i = 0; i < rhs.size(); ++i)
-      rhs[i] = rho_values[i];
+      for (size_t i = 0; i < points.size(); ++i)
+        x[i] = points[i][0];
 
-    poisson.set_rhs(rhs);
+      return eval_rho(it, x, grid_versions, phi_history, Parameters::NV);
+    });
 
     poisson.solve_step();