diff --git a/libnufi_lib.a b/libnufi_lib.a
index fb38e5b..4dbcb62 100644
Binary files a/libnufi_lib.a and b/libnufi_lib.a differ
diff --git a/nufi/fields.h b/nufi/fields.h
index 97aa7cf..0e5b1fe 100644
--- a/nufi/fields.h
+++ b/nufi/fields.h
@@ -76,8 +76,8 @@ inline double eval(double x, const PoissonProblem<1> &poisson,
 
   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));
+  return eval_point_grad<1>(poisson.get_mapping(), poisson.get_dof_handler(),
+                            solution, Point<1>(x));
 }
 
 inline double integral_space_vector(const PoissonProblem<1> &poisson,
diff --git a/nufi/parameters.h b/nufi/parameters.h
index 3b533f8..c379718 100644
--- a/nufi/parameters.h
+++ b/nufi/parameters.h
@@ -36,7 +36,7 @@ constexpr double DT = 1. / 8.;
 constexpr unsigned int TMAX = 100;
 
 // Plotting options
-constexpr int PLOT_FREQUENCY = 10;
+constexpr int PLOT_FREQUENCY = 4;
 constexpr size_t PLOT_NX = CALC_NX;
 constexpr double PLOT_DX = LX / PLOT_NX;
 } // namespace Parameters
diff --git a/nufi/poisson_problem.h b/nufi/poisson_problem.h
index 573ebf2..354903e 100644
--- a/nufi/poisson_problem.h
+++ b/nufi/poisson_problem.h
@@ -33,11 +33,14 @@
 #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 <deal.II/numerics/vector_tools_point_gradient.h>
+#include <deal.II/numerics/vector_tools_point_value.h>
 #include <functional>
 #include <memory>
 #include <string>
@@ -170,14 +173,26 @@ 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);
+double
+eval_point_grad(const Mapping<dim> &mapping, const DoFHandler<dim> &dof_handler,
+                const Vector<double> &solution, const Point<dim> &point) {
+
+  Tensor<1, dim> grad =
+      VectorTools::point_gradient<dim>(mapping, dof_handler, solution, point);
+  double Ex = grad[0];
+
+  return Ex;
 }
 
-// dealii Poisson
+template <int dim>
+double eval_point_value(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);
+}
+// dealii Poisson
 template <int dim> void PoissonProblem<dim>::create_mesh() {
 
   GridGenerator::hyper_cube(triangulation, Parameters::X_DOMAIN_LEFT,
diff --git a/src/main.cc b/src/main.cc
index 66caca3..e5f3c67 100644
--- a/src/main.cc
+++ b/src/main.cc
@@ -1,7 +1,7 @@
 #include <filesystem>
 #include <iostream>
-
 #include <nufi/nufi_solver.h>
+#include <omp.h>
 
 void clear_results_directory(const std::string &dir) {
   if (!std::filesystem::exists(dir) || !std::filesystem::is_directory(dir))
@@ -16,13 +16,13 @@ void clear_results_directory(const std::string &dir) {
 }
 
 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;
diff --git a/src/nufi_solver.cc b/src/nufi_solver.cc
index ac90569..d776d61 100644
--- a/src/nufi_solver.cc
+++ b/src/nufi_solver.cc
@@ -80,7 +80,7 @@ double NuFISolver::eval_rho(unsigned int n, const double x,
                             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;
+  const double v_min = Parameters::V_DOMAIN_LEFT + 0.5 * dv;
 
   double integral = 0.0;
 
@@ -146,14 +146,16 @@ void NuFISolver::run() {
 
     phi_history.push_back(poisson.get_solution());
 
-    std::vector<double> sampled_potential =
-        poisson.sample_electric_potential(x_min, x_max, Nx); // Solution of FE
+    // std::vector<double> sampled_potential =
+    //     poisson.sample_electric_potential(x_min, x_max, Nx); // Solution of
+    //     FE
 
     double timer_elapsed = timer.elapsed();
     double step_time = timer_elapsed - time_elapsed_before;
     total_time += timer_elapsed;
 
     time_file << it << " " << step_time << " " << total_time << "\n";
+    time_file.flush();
 
     std::cout << "step made in " << step_time << " seconds\n\n";
     if (it % Parameters::PLOT_FREQUENCY == 0) {