diff --git a/fields.hpp b/fields.hpp
index 4d8ee0b..fa78be3 100644
--- a/fields.hpp
+++ b/fields.hpp
@@ -20,6 +20,7 @@ inline double f0(const double x,
   return prefactor * gaussian;
 }
 
+
 inline double compute_rho(const double x,
                           const unsigned int Nv = Parameters::NV)
 {
@@ -37,7 +38,7 @@ inline double compute_rho(const double x,
 }
 
 template <int dim>
-class ChargeDensity : public Function<dim>
+class ChargeDensity : public Function<dim> // only uses f0
 {
 public:
   ChargeDensity(double eps,
diff --git a/nufi_poisson.cc b/nufi_poisson.cc
index 2d1dfae..c13ab22 100644
--- a/nufi_poisson.cc
+++ b/nufi_poisson.cc
@@ -3,15 +3,21 @@
 #include "parameters.hpp"
 #include "poisson_problem.hpp"
 
+#include "nufi_solver.hpp"
+
 int main()
 {
   try
   {
-    PoissonProblem<Parameters::DIMENSION> poisson_problem(Parameters::FE_DEGREE,
-                                                          Parameters::NV);
-
-    poisson_problem.run();
+    // PoissonProblem<Parameters::DIMENSION> poisson_problem(Parameters::FE_DEGREE,
+    //                                                       Parameters::NV);
+    //
+    // poisson_problem.run();
+    
+    NuFISolver solver;
+    solver.run();
   }
+
   catch (std::exception &exc)
   {
     std::cerr << std::endl
diff --git a/nufi_solver.hpp b/nufi_solver.hpp
new file mode 100644
index 0000000..68017c6
--- /dev/null
+++ b/nufi_solver.hpp
@@ -0,0 +1,133 @@
+#ifndef NUFI_SOLVER_HPP
+#define NUFI_SOLVER_HPP
+
+#include <deal.II/base/point.h>
+#include <deal.II/base/tensor.h>
+#include <deal.II/numerics/fe_field_function.h>
+
+#include <vector>
+#include <cstddef>
+
+#include "parameters.hpp"
+#include "poisson_problem.hpp"
+#include "fields.hpp" // holds f0(x,v), and compute_rho(x)
+
+using namespace dealii;
+
+class NuFISolver
+{
+public:
+  NuFISolver();
+
+  void run();
+  double eval_rho(size_t n, double x, double u);
+
+private:
+
+  void compute_density();
+  void solve_poisson();
+  void update_distribution();
+
+  double evaluate_E(double x);
+
+  PoissonProblem<1> poisson;
+
+  std::vector<double> coeffs;
+  std::vector<double> rho;
+
+  unsigned int Nt;
+  unsigned int Nx;
+
+  double Lx = Parameters::LX;
+
+  unsigned int order;
+
+  double dt;
+  size_t stride_t;
+};
+
+inline double NuFISolver::evaluate_E(double x)
+{
+  // Wrap x into the periodic domain
+  double x_periodic = x - Lx * std::floor(x / Lx);
+  Point<1> p(x_periodic);
+
+  Functions::FEFieldFunction<1> E_field(
+        poisson.get_dof_handler(),
+        poisson.get_solution()
+    );
+
+  double E_val = 0.0;
+
+  try
+  {
+      // Evaluate the electric field at point p
+      // If your solution represents phi, take negative gradient
+      Tensor<1,1> grad = E_field.gradient(p);
+      E_val = -grad[0];  // -∂φ/∂x
+  }
+  catch (const VectorTools::ExcPointNotAvailableHere &)
+  {
+      // This happens if p lies in an artificial cell in parallel
+      AssertThrow(false, ExcMessage("Point not available on this process."));
+  }
+
+  return E_val;
+}
+
+inline double NuFISolver::eval_rho(size_t n,
+                                      double x,
+                                      double u)
+{
+    if (n == 0)
+        return f0(x,u);
+   
+    double Ex;
+
+    // Initial half-step.
+    Ex = evaluate_E(x);
+    u += 0.5*dt*Ex;
+
+    while ( --n )
+    {
+        x -= dt*u;
+        Ex = evaluate_E(x);
+        u += dt*Ex;
+    }
+
+    // Final half-step.
+    x -= dt*u;
+    Ex = evaluate_E(x);
+    u += 0.5*dt*Ex; // is this line useless ?
+    
+    double x_periodic = x - Lx * std::floor(x / Lx);
+
+    return compute_rho(x_periodic); // compute_rho (from fields.hpp) uses f0
+}
+
+inline void NuFISolver::run()
+{
+    rho.resize(Nx, 0.0);               // initialize density array
+
+    // Main time-stepping loop
+    for (size_t n = 0; n < Nt; ++n)
+    {
+
+        // Solve this logic! To use on deal.ii  
+
+        // 1. Compute charge density rho from current distribution
+        compute_density();
+
+        // 2. Solve Poisson's equation to update electric field
+        solve_poisson();
+
+        // 3. Update distribution function along characteristics
+        update_distribution();
+
+        // Optional: compute ftilda at current step if needed
+        // for demonstration: evaluate ftilda at midpoint x, u = 0
+        // double ft = eval_ftilda(n, 0.5 * poisson.get_Lx(), 0.0);
+    }
+}
+
+#endif
diff --git a/parameters.hpp b/parameters.hpp
index b7c7add..fd45e75 100644
--- a/parameters.hpp
+++ b/parameters.hpp
@@ -1,12 +1,14 @@
 #ifndef PARAMETERS_HPP
 #define PARAMETERS_HPP
 
+#include <cstdlib>
 namespace Parameters
 {
   constexpr unsigned int DIMENSION = 1;
 
   constexpr double X_DOMAIN_LEFT  = 0.0;
   constexpr double X_DOMAIN_RIGHT = 12.0;
+  constexpr double LX = std::abs(X_DOMAIN_RIGHT- X_DOMAIN_LEFT);
 
   constexpr double V_DOMAIN_LEFT = -6.0;
   constexpr double V_DOMAIN_RIGHT = 6.0;
diff --git a/poisson_problem.hpp b/poisson_problem.hpp
index 9c05ab7..086ba07 100644
--- a/poisson_problem.hpp
+++ b/poisson_problem.hpp
@@ -2,7 +2,6 @@
 #define POISSON_PROBLEM_HPP
 
 #include <fstream>
-#include <iostream>
 
 #include <deal.II/base/quadrature_lib.h>
 #include <deal.II/base/logstream.h>
@@ -48,6 +47,9 @@ public:
 
   void set_Nv(unsigned int new_Nv);
 
+  const Vector<double> &get_solution() const { return solution; }
+  const DoFHandler<dim> &get_dof_handler() const { return dof_handler; }
+
 private:
   void create_mesh();
   void setup_system();
@@ -288,6 +290,7 @@ void PoissonProblem<dim>::output_results() const
 template <int dim>
 void PoissonProblem<dim>::run()
 {
+  set_Nv(Parameters::NV); // dont use anywhere else! Other functions still use Parameters::NV.
   create_mesh();
   setup_system();
   assemble_system();