diff --git a/README.md b/README.md
index 74c5412..ff1ba00 100644
--- a/README.md
+++ b/README.md
@@ -3,8 +3,3 @@
 This simulation of the Vlasov-Poisson system in 1x1v dimensions uses
 - [NuFI algorithm](https://doi.org/10.1002/pamm.202300162)
 - [deal.ii](https://dealii.org/) FEM package
-
---- 
-Todo:
-- fix eval and saving fields.
-    - something with periodicity or eval range in x
diff --git a/nufi/fields.h b/nufi/fields.h
index 0e5b1fe..252e819 100644
--- a/nufi/fields.h
+++ b/nufi/fields.h
@@ -4,8 +4,10 @@
 #include "nufi/parameters.h"
 #include "nufi/poisson_problem.h"
 #include <cmath>
+#include <cstddef>
 #include <deal.II/base/function.h>
 #include <deal.II/base/point.h>
+#include <vector>
 
 using namespace dealii;
 
@@ -69,15 +71,22 @@ inline double f0(const double x, const double v,
 }
 
 // wrapper for eval_point() { VectorTools::point_values() }
-inline double eval(double x, const PoissonProblem<1> &poisson,
-                   const Vector<double> &solution) noexcept {
+inline std::vector<double> eval(std::vector<double> &X,
+                                const PoissonProblem<1> &poisson,
+                                const Vector<double> &solution) noexcept {
+  size_t x_size = X.size();
+  std::vector<double> evals(x_size);
+  std::vector<Point<1>> Points(x_size);
 
-  x -= Parameters::X_DOMAIN_LEFT;
+  for (size_t i = 0; i < x_size; ++i) {
+    X[i] = X[i] - Parameters::X_DOMAIN_LEFT;
+    X[i] = X[i] - Parameters::LX * std::floor(X[i] * Parameters::LX_INV);
 
-  x = x - Parameters::LX * std::floor(x * Parameters::LX_INV); // in domain
+    Points[i][0] = X[i];
+  }
 
-  return eval_point_grad<1>(poisson.get_mapping(), poisson.get_dof_handler(),
-                            solution, Point<1>(x));
+  return eval_vector_grad(poisson.get_mapping(), poisson.get_dof_handler(),
+                          solution, Points);
 }
 
 inline double integral_space_vector(const PoissonProblem<1> &poisson,
diff --git a/nufi/nufi_solver.h b/nufi/nufi_solver.h
index 95fc301..5dc46fa 100644
--- a/nufi/nufi_solver.h
+++ b/nufi/nufi_solver.h
@@ -19,16 +19,18 @@ public:
   NuFISolver();
 
   void run();
-  double eval_rho(unsigned int n, const double x,
-                  const PoissonProblem<1> &poisson,
-                  const std::vector<Vector<double>> &phi_history,
-                  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;
-  double eval_f(unsigned int n, double x, double u,
-                const PoissonProblem<1> &poisson,
-                const std::vector<Vector<double>> &phi_history) const;
+  std::vector<double> eval_rho(unsigned int n, std::vector<double> &x,
+                               const PoissonProblem<1> &poisson,
+                               const std::vector<Vector<double>> &phi_history,
+                               const unsigned int Nv = Parameters::NV) const;
+  std::vector<double>
+  eval_ftilda(unsigned int, std::vector<double> &x, double u,
+              const PoissonProblem<1> &poisson,
+              const std::vector<Vector<double>> &phi_history) const;
+  std::vector<double>
+  eval_f(unsigned int n, std::vector<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);
diff --git a/nufi/poisson_problem.h b/nufi/poisson_problem.h
index 354903e..fecbe53 100644
--- a/nufi/poisson_problem.h
+++ b/nufi/poisson_problem.h
@@ -184,6 +184,20 @@ eval_point_grad(const Mapping<dim> &mapping, const DoFHandler<dim> &dof_handler,
   return Ex;
 }
 
+template <int dim>
+std::vector<double> eval_vector_grad(const Mapping<dim> &mapping,
+                                     const DoFHandler<dim> &dof_handler,
+                                     const Vector<double> &solution,
+                                     const std::vector<Point<dim>> &points) {
+  size_t p_size = points.size();
+
+  std::vector<double> Ex(p_size);
+  for (size_t i = 0; i < p_size; ++i)
+    Ex[i] = eval_point_grad(mapping, dof_handler, solution, points[i]);
+
+  return Ex;
+}
+
 template <int dim>
 double eval_point_value(const Mapping<dim> &mapping,
                         const DoFHandler<dim> &dof_handler,
diff --git a/src/nufi_solver.cc b/src/nufi_solver.cc
index d776d61..7a512be 100644
--- a/src/nufi_solver.cc
+++ b/src/nufi_solver.cc
@@ -22,72 +22,127 @@
 
 using namespace dealii;
 
-double
-NuFISolver::eval_ftilda(unsigned int n, double x, double u,
+std::vector<double>
+NuFISolver::eval_ftilda(unsigned int n, std::vector<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;
+  size_t x_size = X.size();
+
+  std::vector<double> U(x_size, u);
+  std::vector<double> results(x_size);
+  if (n == 0) {
+    for (size_t i = 0; i < x_size; ++i)
+      results[i] = f0(X[i], U[i]);
+    return results;
+  }
+
+  std::vector<double> Ex(x_size);
+  std::vector<double> tmp(x_size);
 
   // We omit the initial half-step.
-
   while (--n) {
-    x = x - Parameters::DT * u;
-    Ex = -eval(x, poisson, phi_history[n]);
-    u = u + Parameters::DT * Ex;
+    for (size_t i = 0; i < x_size; ++i)
+      X[i] = X[i] - Parameters::DT * U[i];
+
+    tmp = eval(X, poisson, phi_history[n]); // call eval only once
+
+    for (size_t i = 0; i < x_size; ++i) {
+      Ex[i] = -tmp[i];
+      U[i] = U[i] + Parameters::DT * Ex[i];
+    }
   }
 
   // The final half-step.
-  x = x - Parameters::DT * u;
-  Ex = -eval(x, poisson, phi_history[n]);
-  u += 0.5 * Parameters::DT * Ex;
+  for (size_t i = 0; i < x_size; ++i)
+    X[i] = X[i] - Parameters::DT * U[i];
 
-  return f0(x, u);
+  tmp = eval(X, poisson, phi_history[n]); // call eval only once
+
+  for (size_t i = 0; i < x_size; ++i) {
+    Ex[i] = -tmp[i];
+    U[i] = U[i] + 0.5 * Parameters::DT * Ex[i];
+  }
+  for (size_t i = 0; i < x_size; ++i)
+    results[i] = f0(X[i], U[i]);
+  return results;
 }
 
-double
-NuFISolver::eval_f(unsigned int n, double x, double u,
+std::vector<double>
+NuFISolver::eval_f(unsigned int n, std::vector<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;
+  size_t x_size = X.size();
 
+  std::vector<double> U(x_size, u);
+  std::vector<double> results(x_size);
+  if (n == 0) {
+    for (size_t i = 0; i < x_size; ++i)
+      results[i] = f0(X[i], U[i]);
+    return results;
+  }
+
+  std::vector<double> Ex(x_size);
+
+  std::vector<double> tmp(x_size);
   // Initial half-step.
-  Ex = -eval(x, poisson, phi_history[n]);
-  u += 0.5 * Parameters::DT * Ex;
+  tmp = eval(X, poisson, phi_history[n]); // call eval only once
+  for (size_t i = 0; i < x_size; ++i) {
+    Ex[i] = -tmp[i];
+    U[i] = U[i] + 0.5 * Parameters::DT * Ex[i];
+  }
 
   while (--n) {
-    x = x - Parameters::DT * u;
-    Ex = -eval(x, poisson, phi_history[n]);
-    u = u + Parameters::DT * Ex;
+    for (size_t i = 0; i < x_size; ++i)
+      X[i] = X[i] - Parameters::DT * U[i];
+
+    tmp = eval(X, poisson, phi_history[n]); // call eval only once
+    for (size_t i = 0; i < x_size; ++i) {
+      Ex[i] = -tmp[i];
+      U[i] = U[i] + Parameters::DT * Ex[i];
+    }
   }
 
   // The final half-step.
-  x = x - Parameters::DT * u;
-  Ex = -eval(x, poisson, phi_history[n]);
-  u += 0.5 * Parameters::DT * Ex;
+  for (size_t i = 0; i < x_size; ++i)
+    X[i] = X[i] - Parameters::DT * U[i];
 
-  return f0(x, u);
+  tmp = eval(X, poisson, phi_history[n]); // call eval only once
+  for (size_t i = 0; i < x_size; ++i) {
+    Ex[i] = -tmp[i];
+    U[i] = U[i] + 0.5 * Parameters::DT * Ex[i];
+  }
+
+  for (size_t i = 0; i < x_size; ++i)
+    results[i] = f0(X[i], U[i]);
+  return results;
 }
 
-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 {
+std::vector<double>
+NuFISolver::eval_rho(unsigned int n, std::vector<double> &X,
+                     const PoissonProblem<1> &poisson,
+                     const std::vector<Vector<double>> &phi_history,
+                     const unsigned int Nv) const {
+  size_t x_size = X.size();
   const double dv =
       (Parameters::V_DOMAIN_RIGHT - Parameters::V_DOMAIN_LEFT) / Nv;
   const double v_min = Parameters::V_DOMAIN_LEFT + 0.5 * dv;
 
-  double integral = 0.0;
+  std::vector<double> integral(x_size, 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, phi_history);
-  return 1.0 - integral * dv;
+  std::vector<double> tmp_int(x_size);
+
+  for (unsigned int i = 0; i < Nv; ++i) {
+    tmp_int = eval_ftilda(n, X, v_min + i * dv, poisson,
+                          phi_history); // used eval_ftilda once per i
+    for (size_t ii = 0; ii < x_size; ++ii)
+      integral[ii] += tmp_int[ii];
+  }
+
+  for (size_t i = 0; i < x_size; ++i)
+    integral[i] = 1 - integral[i] * dv;
+  return integral;
 }
 
 void NuFISolver::run() {