diff --git a/fields.hpp b/fields.hpp
index fa78be3..1e15714 100644
--- a/fields.hpp
+++ b/fields.hpp
@@ -33,7 +33,6 @@ inline double compute_rho(const double x,
     const double v = Parameters::V_DOMAIN_LEFT + (i + 0.5) * dv;
     integral += f0(x, v) * dv;
   }
-
   return 1.0 - integral;
 }
 
diff --git a/nufi_poisson.cc b/nufi_poisson.cc
index c13ab22..44e87bf 100644
--- a/nufi_poisson.cc
+++ b/nufi_poisson.cc
@@ -1,8 +1,4 @@
 #include <iostream>
-
-#include "parameters.hpp"
-#include "poisson_problem.hpp"
-
 #include "nufi_solver.hpp"
 
 int main()
diff --git a/nufi_solver.hpp b/nufi_solver.hpp
index 68017c6..3791510 100644
--- a/nufi_solver.hpp
+++ b/nufi_solver.hpp
@@ -1,6 +1,13 @@
+/*
+Todo: 
+- update Nx between timesteps to account for adaptivity changes because Vector rho needs to be resized
+ */
+
 #ifndef NUFI_SOLVER_HPP
 #define NUFI_SOLVER_HPP
 
+#include <cmath>
+#include <cstdlib>
 #include <deal.II/base/point.h>
 #include <deal.II/base/tensor.h>
 #include <deal.II/numerics/fe_field_function.h>
@@ -20,30 +27,28 @@ public:
   NuFISolver();
 
   void run();
-  double eval_rho(size_t n, double x, double u);
+  double eval_rho(unsigned int n, double x, unsigned int Nv = Parameters::NV);
 
 private:
 
-  void compute_density();
-  void solve_poisson();
-  void update_distribution();
+  double eval_ftilda(unsigned int n, double x, double u);
+  void solve_poisson(unsigned int n);
 
   double evaluate_E(double x);
-
-  PoissonProblem<1> poisson;
-
-  std::vector<double> coeffs;
+  
   std::vector<double> rho;
 
-  unsigned int Nt;
+  unsigned int Nt = std::floor(Parameters::TMAX/Parameters::DT);
   unsigned int Nx;
 
   double Lx = Parameters::LX;
 
   unsigned int order;
 
-  double dt;
-  size_t stride_t;
+  double dt = Parameters::DT;
+
+  PoissonProblem<1> poisson;
+
 };
 
 inline double NuFISolver::evaluate_E(double x)
@@ -75,59 +80,117 @@ inline double NuFISolver::evaluate_E(double x)
   return E_val;
 }
 
-inline double NuFISolver::eval_rho(size_t n,
+inline double NuFISolver::eval_ftilda(unsigned int n,
                                       double x,
                                       double u)
 {
-    if (n == 0)
-        return f0(x,u);
-   
-    double Ex;
+  double Lu = std::abs(Parameters::V_DOMAIN_LEFT - Parameters::V_DOMAIN_RIGHT);
 
-    // Initial half-step.
-    Ex = evaluate_E(x);
-    u += 0.5*dt*Ex;
+  if (n == 0)
+      return f0(x, u);
+ 
+  double Ex;
 
-    while ( --n )
-    {
-        x -= dt*u;
-        Ex = evaluate_E(x);
-        u += dt*Ex;
-    }
+  // Initial half-step.
+  Ex = evaluate_E(x);
+  u += 0.5*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);
+  while ( --n )
+  {
+      x -= dt*u;
+      Ex = evaluate_E(x);
+      u += dt*Ex;
+  }
 
-    return compute_rho(x_periodic); // compute_rho (from fields.hpp) uses f0
+  // 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);
+  double u_periodic = u - Lu * std::floor(u / Lu);
+
+  return f0(x_periodic, u_periodic);
+}
+
+inline double NuFISolver::eval_rho(const unsigned int n,
+                          const double x,
+                          const unsigned int Nv)
+{
+  const double dv = (Parameters::V_DOMAIN_RIGHT - Parameters::V_DOMAIN_LEFT) / Nv;
+
+  double integral = 0.0;
+
+  for (unsigned int i = 0; i < Nv; ++i)
+  {
+    const double v = Parameters::V_DOMAIN_LEFT + (i + 0.5) * dv;
+    integral += eval_ftilda(n, x, v) * dv;
+  }
+
+  return 1.0 - integral;
+}
+
+class ChargeDensity_NuFI : public Function<1>
+{
+public:
+  ChargeDensity_NuFI(NuFISolver &solver, size_t n)
+      : solver(solver), n(n) {}
+
+  virtual double value(const Point<1> &p,
+                       [[maybe_unused]] const unsigned int component = 0) const override
+  {
+      double x = p[0];
+
+      return solver.eval_rho(n, x); 
+      // u not used anymore
+  }
+
+private:
+  NuFISolver &solver;
+  size_t n;
+};
+
+
+inline void NuFISolver::solve_poisson(unsigned int n)
+{
+  ChargeDensity_NuFI rho_function(*this, n);
+
+  poisson.set_rhs_function(rho_function);
+
+  poisson.solve_step();
 }
 
 inline void NuFISolver::run()
 {
-    rho.resize(Nx, 0.0);               // initialize density array
+  std::cout << "Starting NuFI solver\n";
 
-    // Main time-stepping loop
-    for (size_t n = 0; n < Nt; ++n)
+  for (unsigned int n = 0; n < Nt; ++n)
+  {
+    std::cout << "Timestep " << n << " / " << Nt << std::endl;
+
+
+    double dx = Lx / Nx;
+
+    for (unsigned int i = 0; i < Nx; ++i)
     {
-
-        // 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);
+      double x = (i + 0.5) * dx;
+      rho[i] = eval_rho(n, x);
     }
+
+    solve_poisson(n);
+  }
+
+  std::cout << "NuFI simulation finished.\n";
 }
 
+inline NuFISolver::NuFISolver()
+  : order(Parameters::FE_DEGREE),
+    poisson(order, Parameters::NV)
+{
+  poisson.initialize();
+
+  Nx = poisson.get_dof_handler().n_dofs();
+  rho.resize(Nx, 0.0);
+
+}
 #endif
diff --git a/parameters.hpp b/parameters.hpp
index fd45e75..97d9fe9 100644
--- a/parameters.hpp
+++ b/parameters.hpp
@@ -20,6 +20,9 @@ namespace Parameters
 
   constexpr double EPS = 0.01;
   constexpr double WAVE_NR = 0.5;
+
+  constexpr double DT=0.05;
+  constexpr unsigned int TMAX = 10;
 }
 
 #endif
diff --git a/poisson_problem.hpp b/poisson_problem.hpp
index 086ba07..b5b6f57 100644
--- a/poisson_problem.hpp
+++ b/poisson_problem.hpp
@@ -1,6 +1,7 @@
 #ifndef POISSON_PROBLEM_HPP
 #define POISSON_PROBLEM_HPP
 
+#include <deal.II/base/function.h>
 #include <fstream>
 
 #include <deal.II/base/quadrature_lib.h>
@@ -43,9 +44,13 @@ class PoissonProblem
 {
 public:
   PoissonProblem(unsigned int degree, unsigned int Nv);
+
+  void initialize();
+  void solve_step();
   void run();
 
   void set_Nv(unsigned int new_Nv);
+  void set_rhs_function(const Function<dim> &rhs);
 
   const Vector<double> &get_solution() const { return solution; }
   const DoFHandler<dim> &get_dof_handler() const { return dof_handler; }
@@ -69,6 +74,8 @@ private:
   Vector<double> solution;      // phi
   Vector<double> system_rhs;
 
+  const Function<dim> *rhs_function;
+
   MappingQ<dim> mapping;
 
   unsigned int Nv;
@@ -81,6 +88,12 @@ void PoissonProblem<dim>::set_Nv(unsigned int new_Nv)
   Nv = new_Nv;
 }
 
+template <int dim>
+void PoissonProblem<dim>::set_rhs_function(const Function<dim> &rhs)
+{
+  rhs_function = &rhs;
+}
+
 template <int dim>
 PoissonProblem<dim>::PoissonProblem(unsigned int degree, unsigned int Nv)
   : fe(degree)
@@ -185,7 +198,7 @@ void PoissonProblem<dim>::assemble_system()
   Vector<double>     cell_rhs(dofs_per_cell);
   std::vector<types::global_dof_index> local_dof_indices(dofs_per_cell);
 
-  ChargeDensity<dim> rhs_function(Parameters::EPS, Parameters::WAVE_NR, Nv);
+  Assert(rhs_function != nullptr, ExcMessage("RHS function not set"));
 
   for (const auto &cell : dof_handler.active_cell_iterators())
   {
@@ -195,7 +208,7 @@ void PoissonProblem<dim>::assemble_system()
 
     for (unsigned int q = 0; q < n_q_points; ++q)
     {
-      const double rho = rhs_function.value(fe_values.quadrature_point(q));
+      const double rho = rhs_function->value(fe_values.quadrature_point(q));
 
       for (unsigned int i = 0; i < dofs_per_cell; ++i)
       {
@@ -286,7 +299,27 @@ void PoissonProblem<dim>::output_results() const
   data_out_E.write_vtk(out2);
 }
 
+template <int dim>
+void PoissonProblem<dim>::initialize()
+{
+  set_Nv(Parameters::NV);
 
+  create_mesh();      // build grid
+  setup_system();     // distribute DoFs and matrices
+}
+
+template <int dim>
+void PoissonProblem<dim>::solve_step()
+{
+  system_matrix = 0;
+  system_rhs = 0;
+
+  assemble_system();
+  solve();
+}
+
+
+// NuFI doesnt use this, kept only for testing. 
 template <int dim>
 void PoissonProblem<dim>::run()
 {