diff --git a/CMakeLists.txt b/CMakeLists.txt
index 54666d7..a565c27 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -23,7 +23,7 @@ find_package(OpenMP REQUIRED)
 add_library(nufi_lib
   src/nufi_solver.cc
   src/save_results.cc
-  src/blas.cc
+  #  src/blas.cc
 )
 
 target_include_directories(nufi_lib PUBLIC
diff --git a/Makefile b/Makefile
index 1b6e392..38f5543 100644
--- a/Makefile
+++ b/Makefile
@@ -142,30 +142,6 @@ nufi_poisson/fast:
 	$(MAKE) $(MAKESILENT) -f CMakeFiles/nufi_poisson.dir/build.make CMakeFiles/nufi_poisson.dir/build
 .PHONY : nufi_poisson/fast
 
-src/blas.o: src/blas.cc.o
-.PHONY : src/blas.o
-
-# target to build an object file
-src/blas.cc.o:
-	$(MAKE) $(MAKESILENT) -f CMakeFiles/nufi_lib.dir/build.make CMakeFiles/nufi_lib.dir/src/blas.cc.o
-.PHONY : src/blas.cc.o
-
-src/blas.i: src/blas.cc.i
-.PHONY : src/blas.i
-
-# target to preprocess a source file
-src/blas.cc.i:
-	$(MAKE) $(MAKESILENT) -f CMakeFiles/nufi_lib.dir/build.make CMakeFiles/nufi_lib.dir/src/blas.cc.i
-.PHONY : src/blas.cc.i
-
-src/blas.s: src/blas.cc.s
-.PHONY : src/blas.s
-
-# target to generate assembly for a file
-src/blas.cc.s:
-	$(MAKE) $(MAKESILENT) -f CMakeFiles/nufi_lib.dir/build.make CMakeFiles/nufi_lib.dir/src/blas.cc.s
-.PHONY : src/blas.cc.s
-
 src/main.o: src/main.cc.o
 .PHONY : src/main.o
 
@@ -248,9 +224,6 @@ help:
 	@echo "... rebuild_cache"
 	@echo "... nufi_lib"
 	@echo "... nufi_poisson"
-	@echo "... src/blas.o"
-	@echo "... src/blas.i"
-	@echo "... src/blas.s"
 	@echo "... src/main.o"
 	@echo "... src/main.i"
 	@echo "... src/main.s"
diff --git a/libnufi_lib.a b/libnufi_lib.a
index 14ea604..f66070d 100644
Binary files a/libnufi_lib.a and b/libnufi_lib.a differ
diff --git a/nufi/blas.h b/nufi/blas.h
deleted file mode 100644
index 2f90e84..0000000
--- a/nufi/blas.h
+++ /dev/null
@@ -1,54 +0,0 @@
-#ifndef NUFI_BLAS_H
-#define NUFI_BLAS_H
-
-#include <cstddef>
-
-/*!
- * \brief Convenience wrappers for BLAS, with overloads for single and double
- *        precision.
- */
-namespace blas
-{
-
-double dot( const size_t n, const double *x, size_t incx,
-                            const double *y, size_t incy );
-
-float  dot( const size_t n, const float  *x, size_t incx,
-                            const float  *y, size_t incy );
-
-void axpy( size_t n, double alpha, const double *x, size_t incx,
-                                         double *y, size_t incy );
-
-void axpy( size_t n, float  alpha, const float  *x, size_t incx,
-                                         float  *y, size_t incy );
-
-
-void scal( size_t n, double alpha, double *x, size_t incx );
-void scal( size_t n, float  alpha, float  *x, size_t incx );
-
-void copy( size_t n, const double *x, size_t incx, double *y, size_t incy );
-void copy( size_t n, const float  *x, size_t incx, float  *y, size_t incy );
-
-void ger( const size_t M, const size_t N, const double alpha,
-          const double *X, const size_t incX, const double *Y, const size_t incY,
-          double *A, const size_t lda);
-
-void ger( const size_t M, const size_t N, const float alpha,
-          const float *X, const size_t incX, const float *Y, const size_t incY,
-          float *A, const size_t lda);
-
-
-void gemv( const char trans, size_t m, size_t n,
-		   double alpha, const double *a, size_t lda,
-           const double *x, size_t incx, double beta,
-           double *y, size_t incy );
-
-void gemv( const char trans, size_t m, size_t n,
-		   float alpha, const float *a, size_t lda,
-           const float *x, size_t incx, float beta,
-           float *y, size_t incy );
-}
-
-
-#endif
-
diff --git a/nufi/fields.h b/nufi/fields.h
index 18e5e58..9268c47 100644
--- a/nufi/fields.h
+++ b/nufi/fields.h
@@ -1,35 +1,79 @@
 #ifndef FIELDS_H
 #define FIELDS_H
 
+#include "nufi/parameters.h"
+#include "poisson_problem.h"
 #include <cmath>
 #include <deal.II/base/function.h>
-#include "nufi/parameters.h"
-#include "nufi/splines.h"
-#include "nufi/lsmr.h"
+#include <deal.II/base/point.h>
 
 using namespace dealii;
 
-inline double f0(const double x,
-                 const double v,
+inline std::vector<int> Indices_of_points(const std::vector<double> &points, double x_min, double x_max, double dx, int grid_type=0)
+{ 
+    // grid type:
+    // 0 => uniform
+    // 1 => non uniform (TODO)
+    
+    if (dx <= 0.0) {
+        throw std::invalid_argument("dx must be positive");
+    }
+    if (x_max <= x_min) {
+        throw std::invalid_argument("x_max must be > x_min");
+    }
+
+    std::vector<int> indices;
+    indices.reserve(points.size());
+
+    switch (grid_type) {
+        case 0:
+            {
+            const double L = x_max - x_min;
+            const int N = std::floor(L/dx);
+
+
+            for (double x : points) //GPT loop, to check
+            {
+                x-= x_min;
+                x = x - L * std::floor(x/L);
+
+                int i = static_cast<int>(std::floor(x / dx));
+
+                // safety: handle rare edge case due to floating precision
+                if (i == N) i = 0;
+
+                indices.push_back(i);
+            }
+            }
+        case 1:
+            {
+            throw std::invalid_argument("Case for non uniform grid is not completed");
+            }
+        default:
+            throw std::invalid_argument("Invalid grid_type argument");
+
+    }
+    return indices;
+}
+
+inline double f0(const double x, const double v,
                  const double eps = Parameters::EPS,
-                 const double k   = Parameters::WAVE_NR)
-{
-  const double prefactor = Parameters::F0_FACTOR * (1.0 + eps * std::cos(k*x));
-  const double gaussian  = v*v  * std::exp(-0.5 * v*v);
+                 const double k = Parameters::WAVE_NR) {
+  const double prefactor =
+      Parameters::F0_FACTOR * (1.0 + eps * std::cos(k * x));
+  const double gaussian = v * v * std::exp(-0.5 * v * v);
 
   return prefactor * gaussian;
 }
 
-
 inline double compute_rho(const double x,
-                          const unsigned int Nv = Parameters::NV)
-{
-  const double dv = (Parameters::V_DOMAIN_RIGHT - Parameters::V_DOMAIN_LEFT) / Nv;
+                          const unsigned int Nv = Parameters::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)
-  {
+  for (unsigned int i = 0; i < Nv; ++i) {
     const double v = Parameters::V_DOMAIN_LEFT + (i + 0.5) * dv;
     integral += f0(x, v) * dv;
   }
@@ -37,185 +81,82 @@ inline double compute_rho(const double x,
   return 1.0 - integral;
 }
 
-
-template <size_t dx = 0>
-double eval(double x, const double *coeffs) noexcept
-{
-    using std::floor;
-
-    // Shift to a box that starts at 0.
-    x -= Parameters::X_DOMAIN_LEFT;
-
-    // Get "periodic position" in box at origin.
-    x = x - Parameters::LX * floor( x*Parameters::LX_INV ); 
-
-    // Knot number
-    double x_knot = floor( x*Parameters::SPLINE_DX_INV); 
-
-    size_t ii = static_cast<size_t>(x_knot);
-
-    // Convert x to reference coordinates.
-    x = x*Parameters::SPLINE_DX_INV - x_knot;
-
-    // Scale according to derivative.
-    double factor = 1;
-    for ( size_t i = 0; i < dx; ++i ) factor *= 1*Parameters::SPLINE_DX_INV;
-
-    return factor*splines1d::eval<double,Parameters::SPLINE_ORDER,dx>(x, coeffs + ii);
+double eval(double x, const PoissonProblem<1> &poisson) noexcept {
+    return poisson.evaluate_potential(Point<1>(x));
 }
 
-template <typename real, size_t order>
-void interpolate( real *coeffs, const real *values)
-{
-    std::unique_ptr<real[]> tmp { new real[ Parameters::SPLINE_NX ] };
-
-    for ( size_t i = 0; i < Parameters::SPLINE_NX; ++i )
-        tmp[ i ] = coeffs[ i ];
-
-    struct mat_t
-    {
-        real  N[ order ];
-
-        mat_t()
-        {
-            splines1d::N<real,order>(0,N);
-        }
-
-        void operator()( const real *in, real *out ) const
-        {
-            #pragma omp parallel for
-            for ( size_t i = 0; i < Parameters::SPLINE_NX; ++i )
-            {
-                real result = 0;
-                if ( i + order <= Parameters::SPLINE_NX )
-                {
-                    for ( size_t ii = 0; ii < order; ++ii )
-                        result += N[ii] * in[ i + ii ];
-                }
-                else
-                {
-                    for ( size_t ii = 0; ii < order; ++ii )
-                        result += N[ii]*in[ (i+ii) % Parameters::SPLINE_NX];
-                }
-                out[ i ] = result;
-            }
-        }
-    };
-
-    struct transposed_mat_t
-    {
-        real  N[ order ];
-
-        transposed_mat_t()
-        {
-            splines1d::N<real,order>(0,N);
-        }
-
-        void operator()( const real *in, real *out ) const
-        {
-            for ( size_t i = 0; i < Parameters::SPLINE_NX; ++i )
-                out[ i ] = 0;
-
-            for ( size_t i = 0; i < Parameters::SPLINE_NX; ++i )
-            {
-                if ( i + order <= Parameters::SPLINE_NX )
-                {
-                    for ( size_t ii = 0; ii < order; ++ii )
-                        out[ i + ii ]  += N[ii] * in[ i ];
-                }
-                else
-                {
-                    for ( size_t ii = 0; ii < order; ++ii )
-                        out[ (i+ii) % Parameters::SPLINE_NX ] += N[ii]*in[ i ];
-                }
-            }
-        }
-    };
-
-    mat_t M; transposed_mat_t Mt;
-    lsmr_options<real> opt; opt.silent = true;
-    lsmr( Parameters::SPLINE_NX, Parameters::SPLINE_NX , M, Mt, values, tmp.get(), opt );
-
-    if ( opt.iter == opt.max_iter )
-        std::cerr << "Warning. LSMR did not converge.\n";
-
-    for ( size_t i = 0; i < Parameters::SPLINE_NX + order - 1; ++i )
-        coeffs[ i ] = tmp[ i % Parameters::SPLINE_NX ];
-}
-
-inline double integral_space_vector(const double *current_coeffs, double dx = Parameters::SPLINE_DX, size_t Nx = Parameters::SPLINE_NX)
-{
+inline double integral_space_vector(const PoissonProblem<1> &poisson,
+                                    double dx = Parameters::PLOT_DX,
+                                    size_t Nx = Parameters::PLOT_NX) {
   double integral = 0.0;
   double xmin = Parameters::X_DOMAIN_LEFT;
-  #pragma omp parallel for reduction (+:integral)
-  for (size_t i=0; i<Nx ; ++i) {
+#pragma omp parallel for reduction(+ : integral)
+  for (size_t i = 0; i < Nx; ++i) {
     double x = xmin + i * dx;
-    integral += eval<1>(x, current_coeffs);
+    integral += eval(x, poisson);
   }
-  return integral*dx;
+  return integral * dx;
 };
 
-inline double integral_space_vector_squared(const double *current_coeffs, double dx = Parameters::SPLINE_DX, size_t Nx = Parameters::SPLINE_NX)
-{
+inline double integral_space_vector_squared(const PoissonProblem<1> &poisson,
+                                            double dx = Parameters::PLOT_DX,
+                                            size_t Nx = Parameters::PLOT_NX) {
   double integral = 0.0;
   double xmin = Parameters::X_DOMAIN_LEFT;
-  #pragma omp parallel for reduction (+:integral)
-  for (size_t i=0; i<Nx ; ++i) {
-    double x = xmin + i*dx;
-    double val = eval<1>(x, current_coeffs); 
-    integral += val*val;
+#pragma omp parallel for reduction(+ : integral)
+  for (size_t i = 0; i < Nx; ++i) {
+    double x = xmin + i * dx;
+    double val = eval(x, poisson);
+    integral += val * val;
   }
-  return integral*dx;
+  return integral * dx;
 };
 
 class Gradient {
 public:
-    Gradient(double xmin, double xmax, unsigned int Nx)
-        : xmin_(xmin), xmax_(xmax), Nx_(Nx)
-    {
-        if (xmax_ <= xmin_) {
-            throw std::invalid_argument("xmax must be greater than xmin");
-        }
+  Gradient(double xmin, double xmax, unsigned int Nx)
+      : xmin_(xmin), xmax_(xmax), Nx_(Nx) {
+    if (xmax_ <= xmin_) {
+      throw std::invalid_argument("xmax must be greater than xmin");
+    }
+  }
+
+  std::vector<double> compute(const std::vector<double> &values) const {
+    size_t n = values.size();
+    if (n < 2) {
+      throw std::invalid_argument("Need at least 2 points");
     }
 
-    std::vector<double> compute(const std::vector<double>& values) const {
-        size_t n = values.size();
-        if (n < 2) {
-            throw std::invalid_argument("Need at least 2 points");
-        }
+    std::vector<double> grad(n);
 
-        std::vector<double> grad(n);
+    double dx = (xmax_ - xmin_) / (n - 1);
+    // periodic boundaries
+    grad[0] = -(values[1] - values[n - 1]) / (2.0 * dx);
+    grad[n - 1] = -(values[0] - values[n - 2]) / (2.0 * dx);
 
-        double dx = (xmax_ - xmin_) / (n-1);
-        // periodic boundaries
-        grad[0]     = -(values[1]   - values[n-1]) / (2.0 * dx);
-        grad[n-1]   = -(values[0]   - values[n-2]) / (2.0 * dx);
-
-        for (size_t i = 1; i < n-1; ++i) {
-            grad[i] = -(values[i+1] - values[i-1]) / (2.0 * dx);
-        }
-
-
-        return grad;
+    for (size_t i = 1; i < n - 1; ++i) {
+      grad[i] = -(values[i + 1] - values[i - 1]) / (2.0 * dx);
     }
+
+    return grad;
+  }
+
 private:
-    double xmin_;
-    double xmax_;
-    [[maybe_unused]] unsigned int Nx_;
+  double xmin_;
+  double xmax_;
+  [[maybe_unused]] unsigned int Nx_;
 };
 
 template <int dim>
 class ChargeDensity : public Function<dim> // only uses f0
 {
 public:
-  ChargeDensity(double eps,
-                double k,
-                unsigned int Nv)
-    : Function<dim>(1), eps(eps), k(k), Nv(Nv) {}
+  ChargeDensity(double eps, double k, unsigned int Nv)
+      : Function<dim>(1), eps(eps), k(k), Nv(Nv) {}
 
-  virtual double value(const Point<dim> &p,
-                       [[maybe_unused]] const unsigned int component = 0) const override
-  {
+  virtual double
+  value(const Point<dim> &p,
+        [[maybe_unused]] const unsigned int component = 0) const override {
     return compute_rho(p[0], Nv);
   }
 
@@ -225,5 +166,4 @@ private:
   const unsigned int Nv;
 };
 
-
 #endif
diff --git a/nufi/lsmr.h b/nufi/lsmr.h
deleted file mode 100644
index 144793f..0000000
--- a/nufi/lsmr.h
+++ /dev/null
@@ -1,252 +0,0 @@
-#ifndef LSMR_H
-#define LSMR_H
-
-#include <cmath>
-#include <limits>
-#include <iomanip>
-#include <iostream>
-#include "nufi/blas.h"
-
-template <typename real>
-struct lsmr_options
-{
-    ///////////
-    // INPUT //
-    ///////////
-
-    // Whether to print messages to std::cout.
-    bool silent            = true;
-
-    // Residual of normal equations AᵀAx = Aᵀb
-    bool relative_residual  = true;
-    real target_residual    = std::numeric_limits<real>::epsilon();
-    size_t      max_iter    = 1000;
-
-    // How many Lánczos vectors to keep for local reorthogonalisation.
-    // Choose zero for no reorthogonalisation, pure LSMR.
-    // Choose a large value for complete reorthognalisation.
-    //
-    // In an ideal world without roundoff errors, this would have no effect
-    // at all, as the Lánczos vectors would be perfectly orthogonal. In practice
-    // this property is lost rather quickly.  One may choose to store some of
-    // the most recent Lánczos vectors to enforce this property manually. This
-    // increase convergence speed at the cost of additional memory requirements.
-    size_t reorthogonalise_u = 50;
-    size_t reorthogonalise_v = 50;
-
-    ////////////
-    // OUTPUT //
-    ////////////
-
-    // Iteration count and reached residual.
-    // Estimates of ‖A‖ and cond(A)
-    size_t iter; real residual;
-    real norm_A_estimate, cond_estimate;
-};
-
-template <typename real, typename mat, typename transposed_mat>
-void lsmr( size_t m, size_t n, const mat& A, const transposed_mat& At,
-           const real *b, real *x, lsmr_options<real> &S );
-
-namespace lsmr_impl
-{
-
-template <typename real>
-real norm( size_t n, const real *x )
-{
-    using std::hypot;
-
-    real result = 0;
-    for ( size_t i = 0; i < n; ++i )
-        result = hypot(result,x[i]);
-
-    return result;
-}
-
-// Reorthognalise u with respect to the previous vectors in buffer,
-// using the modified Gram–Schmidt process. Overwrite the oldest vector
-// in buffer when full.
-template <typename real>
-void reorthogonalise( real *buf, size_t n, size_t buffer_max,  
-                      real *u,   size_t iter )
-{
-    using std::min;
-    using blas::dot;
-    using blas::axpy;
-    using blas::scal;
-    using blas::copy;
-
-    size_t n_buffered = min( iter+1, buffer_max );
-    for ( size_t i = 0; i < n_buffered; ++i )
-    {
-        real fac = -dot( n, u, 1, buf + i*n, 1 );
-        axpy( n, fac, buf + i*n, 1, u, 1 );
-    }
-
-    scal( n, 1/norm(n,u), u, 1 );
-    copy( n, u, 1, buf + ((iter+1)%buffer_max)*n, 1 );
-}
-
-}
-
-template <typename real, typename mat, typename transposed_mat>
-void lsmr( size_t m, size_t n, const mat& A, const transposed_mat& At,
-           const real *b, real *x, lsmr_options<real> &S )
-{
-    using std::min;
-    using std::max;
-    using std::abs;
-    using std::swap;
-    using std::hypot;
-    using blas::axpy; 
-    using blas::scal; 
-    using blas::copy;
-    using lsmr_impl::norm;
-    using lsmr_impl::reorthogonalise;
-
-
-    // Allocation of buffers.
-    size_t max_buf = min(n,m)-1;
-    S.reorthogonalise_u = min(S.reorthogonalise_u,max_buf);
-    S.reorthogonalise_v = min(S.reorthogonalise_v,max_buf);
-    size_t u_buffer_size = max( S.reorthogonalise_u, size_t(1) );
-    size_t v_buffer_size = max( S.reorthogonalise_v, size_t(1) );
-
-    std::unique_ptr<real[]> data { new real[ n*( 4 + v_buffer_size )  +
-                                             m*( 2 + u_buffer_size ) ] {} };
-
-    real *u     = data.get();
-    real *utmp  = u     + m;
-    real *ubuf  = utmp  + m;
-    real *v     = ubuf  + m*u_buffer_size;
-    real *vtmp  = v     + n;
-    real *h     = vtmp  + n;
-    real *h_bar = h     + n;
-    real *vbuf  = h_bar + n;
-
-    At(b,v);
-    const real norm_ATb = norm(n,v);
-
-
-    A(x,u); axpy(m,real(-1),b,1,u,1); 
-    scal(m, real(-1), u, 1 );         // u = b - Ax;
-   
-    real alpha = 0; 
-    real beta  = norm(m,u); 
-
-    if ( beta > real(0) )
-    {
-        scal(m, real(1)/beta, u, 1 ); // u = b - Ax / norm(b-Ax)
-        At(u,v);                      // v = At*u
-        alpha = norm(n,v);
-    }
-
-    if ( alpha > real(0) )
-        scal(n, real(1)/alpha, v, 1 ); // v = At*u/norm(At*u)
-
-    copy(n,u,1,ubuf,1);                // u_buf.col(0) = u_buf
-    copy(n,v,1,vbuf,1);                // v_buf.col(0) = v
-    copy(n,v,1,h,1);                   // h = v
-
-    if ( alpha * beta == real(0) ) return;
-
-    
-    real alpha_bar = alpha, zeta_bar = alpha*beta;
-    real rho = 1, rho_bar = 1, c_bar = 1, s_bar = 0;
-    real c, s, theta, zeta, theta_bar, rho_prev, rho_bar_prev;
-
-    // For estimating the condition number.
-    real   sigma_max = 0,   sigma_min = std::numeric_limits<real>::max();
-    real rho_bar_max = 0, rho_bar_min = std::numeric_limits<real>::max();
-
-    S.norm_A_estimate = 0;
-    for ( S.iter = 0; S.iter < S.max_iter; ++S.iter )
-    {
-        // Continue the bidiagonalisation.
-        A(v,utmp); axpy(m,-alpha,u,1,utmp,1); swap(u,utmp); // u = A*v - alpha*u
-        beta = norm(m,u);
-
-        if ( beta > 0 )
-        {
-            scal(m, real(1)/beta, u, 1 );
-            if ( S.reorthogonalise_u )
-                reorthogonalise( ubuf, m, u_buffer_size, u, S.iter );
-
-            S.norm_A_estimate = hypot( alpha, S.norm_A_estimate );
-            S.norm_A_estimate = hypot( beta , S.norm_A_estimate );
-
-            At(u,vtmp); axpy(n,-beta,v,1,vtmp,1); swap(v,vtmp); // v = At*u - beta*v
-            alpha = norm(n,v);
-
-            if ( alpha > 0 )
-            {
-                scal(n,real(1)/alpha, v, 1 );
-                if ( S.reorthogonalise_v )
-                    reorthogonalise( vbuf, n, v_buffer_size, v, S.iter );
-            }
-        }
-
-        // Construct and apply rotation P_k
-        rho_prev     = rho;
-        rho          = hypot(alpha_bar,beta);
-        c            = alpha_bar/rho;
-        s            = beta/rho;
-        theta        = s*alpha;
-        alpha_bar    = c*alpha;
-
-        // Construct and apply rotation \bar{P}_k
-        rho_bar_prev = rho_bar;
-        if ( S.iter )
-        {
-            rho_bar_max = max( rho_bar, rho_bar_max );
-            rho_bar_min = min( rho_bar, rho_bar_min );
-        }
-        theta_bar    = s_bar*rho;
-        rho_bar      = hypot( c_bar*rho, theta );
-        if ( S.iter )
-        {
-            sigma_max    = max( rho_bar_max, c_bar*rho );
-            sigma_min    = min( rho_bar_min, c_bar*rho );
-        }
-        c_bar        = c_bar * rho/rho_bar;
-        s_bar        = theta/rho_bar;
-        zeta         = c_bar * zeta_bar;
-        zeta_bar     = -s_bar*zeta_bar;
-
-
-        // Update h, h_bar, x
-        scal(n, -(theta_bar*rho)/(rho_prev*rho_bar_prev), h_bar, 1 ) ;
-        axpy(n, real(1), h, 1, h_bar, 1 );             // h_bar = h - factor*h_bar
-
-        axpy( n, zeta/(rho*rho_bar), h_bar, 1, x, 1 ); // x += factor * h_bar
-
-        scal(n, -theta/rho, h, 1 );
-        axpy(n, real(1), v, 1, h, 1 );                 // h = v - factor*h;
-
-        // Estimate quantities.
-        if ( S.relative_residual ) S.residual = abs(zeta_bar)/norm_ATb;
-        else                       S.residual = abs(zeta_bar);
-        S.cond_estimate = sigma_max / sigma_min;
-
-        if ( S.residual <= S.target_residual )
-        {
-            if ( S.silent == false )
-            {
-                std::cout << "LSMR: Iteration:  " << std::setw(4)  << S.iter << ", "
-                          << "Residual:  "        << std::setw(12) << std::scientific << S.residual << ", "
-                          << "cond estimate:  "   << std::setw(12) << std::scientific << S.cond_estimate << ".\n";
-            }
-            return;
-        }
-
-        if ( S.silent == false && (S.iter%10) == 0 )
-        {
-            std::cout << "LSMR: Iteration:  " << std::setw(4)  << S.iter << ", "
-                      << "Residual:  "        << std::setw(12) << std::scientific << S.residual << ", "
-                      << "cond estimate:  "   << std::setw(12) << std::scientific << S.cond_estimate << ".\n";
-        }
-    }
-}
-
-#endif
-
diff --git a/nufi/nufi_solver.h b/nufi/nufi_solver.h
index a023895..522b0b5 100644
--- a/nufi/nufi_solver.h
+++ b/nufi/nufi_solver.h
@@ -20,15 +20,15 @@ public:
   NuFISolver();
 
   void run();
-  double eval_rho(unsigned int n, double x, const double *E_coeffs, unsigned int Nv = Parameters::NV) const;
-  double eval_ftilda(unsigned int n, double x, double u, const double *E_coeffs) const;
-  double eval_f(unsigned int n, double x, double u, const double *E_coeffs) const;
+  double eval_rho(unsigned int n, double x, const PoissonProblem<1> &poisson, unsigned int Nv = Parameters::NV) const;
+  double eval_ftilda(unsigned int n, double x, double u, const PoissonProblem<1> &poisson) const;
+  double eval_f(unsigned int n, double x, double u, const PoissonProblem<1> &poisson) const;
 
 private:
 
 
   unsigned int Nt = std::floor(Parameters::TMAX/Parameters::DT);
-  unsigned int Nx = Parameters::SPLINE_NX;
+  unsigned int Nx = Parameters::CALC_NX;
 
   double Lx = Parameters::LX;
 
diff --git a/nufi/parameters.h b/nufi/parameters.h
index 97618e6..0049d18 100644
--- a/nufi/parameters.h
+++ b/nufi/parameters.h
@@ -13,34 +13,33 @@ namespace Parameters
   constexpr double LX = std::abs(X_DOMAIN_RIGHT- X_DOMAIN_LEFT);
   constexpr double LX_INV = 1/LX;
 
+  constexpr size_t CALC_NX = 256;
+  constexpr double CALC_DX = LX/CALC_NX;
+
   constexpr double V_DOMAIN_LEFT = -10.;
   constexpr double V_DOMAIN_RIGHT = 10.;
 
-  constexpr unsigned int NV = 512;
+  constexpr unsigned int NV = 256;
   constexpr double DV = std::abs(V_DOMAIN_RIGHT - V_DOMAIN_LEFT)/NV;
 
   // deal.ii options
   constexpr unsigned int GLOBAL_REFINEMENT = 8;
-  constexpr unsigned int FE_DEGREE = 4;
-  constexpr unsigned int CONVERGENCE_ITERATIONS = 10000;
-  constexpr double CONVERGENCE_LIMIT = 1e-12;
+  constexpr unsigned int FE_DEGREE = 3;
+  constexpr unsigned int CONVERGENCE_ITERATIONS = 5000;
+  constexpr double CONVERGENCE_LIMIT = 1e-8;
 
   constexpr double EPS = 0.01;
   constexpr double WAVE_NR = 0.5;
   constexpr double F0_FACTOR = 0.39894228040143267793994; // 1/sqrt(2pi)
 
   // NUFI options
-  constexpr double DT=1./16.;
-  constexpr unsigned int TMAX = 500;
-
-  //spline options
-  constexpr int SPLINE_NX = 256;
-  constexpr double SPLINE_DX = LX/(SPLINE_NX);
-  constexpr double SPLINE_DX_INV = 1/SPLINE_DX;
-  constexpr size_t SPLINE_ORDER = 4;
+  constexpr double DT=1./4.;
+  constexpr unsigned int TMAX = 50;
 
   //Plotting options
-  constexpr int PLOT_FREQUENCY = 10;
+  constexpr int PLOT_FREQUENCY = 20;
+  constexpr size_t PLOT_NX = 256;
+  constexpr double PLOT_DX = LX/PLOT_NX;
 }
 
 #endif
diff --git a/nufi/poisson_problem.h b/nufi/poisson_problem.h
index 7571603..3de0c54 100644
--- a/nufi/poisson_problem.h
+++ b/nufi/poisson_problem.h
@@ -3,38 +3,38 @@
 
 #include <deal.II/base/function.h>
 
+#include <deal.II/base/index_set.h>
+#include <deal.II/base/logstream.h>
 #include <deal.II/base/mpi_remote_point_evaluation.h>
 #include <deal.II/base/point.h>
 #include <deal.II/base/quadrature_lib.h>
-#include <deal.II/base/logstream.h>
 #include <deal.II/base/template_constraints.h>
 #include <deal.II/base/tensor.h>
 #include <deal.II/base/utilities.h>
-#include <deal.II/base/index_set.h>
 
-#include <deal.II/lac/vector.h>
-#include <deal.II/lac/full_matrix.h>
-#include <deal.II/lac/sparse_matrix.h>
-#include <deal.II/lac/dynamic_sparsity_pattern.h>
-#include <deal.II/lac/solver_cg.h>
-#include <deal.II/lac/precondition.h>
 #include <deal.II/lac/affine_constraints.h>
+#include <deal.II/lac/dynamic_sparsity_pattern.h>
+#include <deal.II/lac/full_matrix.h>
+#include <deal.II/lac/precondition.h>
+#include <deal.II/lac/solver_cg.h>
+#include <deal.II/lac/sparse_matrix.h>
+#include <deal.II/lac/vector.h>
 
-#include <deal.II/grid/tria.h>
 #include <deal.II/grid/grid_generator.h>
 #include <deal.II/grid/grid_tools.h>
+#include <deal.II/grid/tria.h>
 
 #include <deal.II/dofs/dof_handler.h>
-#include <deal.II/dofs/dof_tools.h>
 #include <deal.II/dofs/dof_renumbering.h>
+#include <deal.II/dofs/dof_tools.h>
 
 #include <deal.II/fe/fe_q.h>
 #include <deal.II/fe/fe_values.h>
 
 #include <deal.II/numerics/data_out.h>
-#include <deal.II/numerics/vector_tools.h>
-#include <deal.II/numerics/matrix_tools.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>
@@ -49,9 +49,7 @@ using namespace dealii;
 
 // =-=-=-=-= Poisson Solver =-=-=-=-=
 
-template <int dim>
-class PoissonProblem
-{
+template <int dim> class PoissonProblem {
 public:
   PoissonProblem(unsigned int degree);
 
@@ -64,8 +62,15 @@ public:
   const Vector<double> &get_solution() const { return solution; }
   const DoFHandler<dim> &get_dof_handler() const { return dof_handler; }
 
-  std::vector<double> sample_electric_field(double x_min, double x_max, unsigned int Nx);
-  std::vector<double> sample_electric_potential(double x_min, double x_max, unsigned int Nx);
+  std::vector<double> sample_electric_field(double x_min, double x_max,
+                                            unsigned int Nx);
+  std::vector<double> sample_electric_potential(double x_min, double x_max,
+                                                unsigned int Nx);
+
+  double evaluate_potential(const Point<dim> &p) const
+{
+  return fe_field_function->value(p);
+}
 
 private:
   void create_mesh();
@@ -74,100 +79,93 @@ private:
   void solve();
 
   Triangulation<dim> triangulation;
-  FE_Q<dim>          fe;
-  DoFHandler<dim>    dof_handler;
+  FE_Q<dim> fe;
+  DoFHandler<dim> dof_handler;
 
   AffineConstraints<double> constraints;
 
-  SparsityPattern      sparsity_pattern;
+  SparsityPattern sparsity_pattern;
   SparseMatrix<double> system_matrix;
 
-  Vector<double> solution;      // phi
+  Vector<double> solution; // phi
   Vector<double> system_rhs;
 
   std::unique_ptr<const Function<dim>> rhs_function;
 
   MappingQ<dim> mapping;
+
+  std::unique_ptr<Functions::FEFieldFunction<dim>> fe_field_function;
 };
 
 // Utilities
 
 template <int dim>
-void PoissonProblem<dim>::set_rhs_function(std::unique_ptr<Function<dim>> rhs)
-{
+void PoissonProblem<dim>::set_rhs_function(std::unique_ptr<Function<dim>> rhs) {
   rhs_function = std::move(rhs);
 }
 
 template <int dim>
 PoissonProblem<dim>::PoissonProblem(unsigned int degree)
-  : fe(degree)
-  , dof_handler(triangulation)
-  , mapping(degree)
-{}
+    : fe(degree), dof_handler(triangulation), mapping(degree) {}
 
 template <int dim>
-std::vector<double> PoissonProblem<dim>::sample_electric_field(double x_min,double x_max,unsigned int Nx)
-{
+std::vector<double>
+PoissonProblem<dim>::sample_electric_field(double x_min, double x_max,
+                                           unsigned int Nx) {
   std::vector<double> E_values(Nx);
 
   const double dx = (x_max - x_min) / (Nx - 1);
 
-  for (unsigned int i = 0; i < Nx; ++i)
-  {
-      const double x = x_min + i * dx;
-      const Point<dim> point(x);
+  for (unsigned int i = 0; i < Nx; ++i) {
+    const double x = x_min + i * dx;
+    const Point<dim> point(x);
 
-      // 1. Find the active cell containing x
-      const auto cell_point_pair =
-          GridTools::find_active_cell_around_point(mapping,
-                                                   dof_handler,
-                                                   point);
+    // 1. Find the active cell containing x
+    const auto cell_point_pair =
+        GridTools::find_active_cell_around_point(mapping, dof_handler, point);
 
-      const auto cell = cell_point_pair.first;
-      const Point<dim> &unit_point = cell_point_pair.second;
+    const auto cell = cell_point_pair.first;
+    const Point<dim> &unit_point = cell_point_pair.second;
 
-      // 2. FEPointEvaluation expects an ArrayView of points
-      std::vector<Point<dim>> points(1, unit_point);
-      ArrayView<const Point<dim>> point_view(points);
+    // 2. FEPointEvaluation expects an ArrayView of points
+    std::vector<Point<dim>> points(1, unit_point);
+    ArrayView<const Point<dim>> point_view(points);
 
-      FEPointEvaluation<1, dim> evaluator(mapping,
-                                          dof_handler.get_fe(),
-                                          update_gradients);
+    FEPointEvaluation<1, dim> evaluator(mapping, dof_handler.get_fe(),
+                                        update_gradients);
 
-      // reinit with ArrayView of points
-      evaluator.reinit(cell, point_view);
+    // reinit with ArrayView of points
+    evaluator.reinit(cell, point_view);
 
-      Vector<double> local_dofs(dof_handler.get_fe().dofs_per_cell);
-      cell->get_dof_values(solution, local_dofs);
+    Vector<double> local_dofs(dof_handler.get_fe().dofs_per_cell);
+    cell->get_dof_values(solution, local_dofs);
 
-      // 3. Evaluate gradient at this point
-      evaluator.evaluate(local_dofs, EvaluationFlags::gradients);
+    // 3. Evaluate gradient at this point
+    evaluator.evaluate(local_dofs, EvaluationFlags::gradients);
 
-      const Tensor<1, dim> grad_phi = evaluator.get_gradient(0);
+    const Tensor<1, dim> grad_phi = evaluator.get_gradient(0);
 
-      // 4. Compute E = -grad(phi)
-      E_values[i] = -grad_phi[0];
+    // 4. Compute E = -grad(phi)
+    E_values[i] = -grad_phi[0];
   }
 
   return E_values;
 }
 
 template <int dim>
-std::vector<double> PoissonProblem<dim>::sample_electric_potential(
-    double x_min,
-    double x_max,
-    unsigned int Nx)
-{
+std::vector<double>
+PoissonProblem<dim>::sample_electric_potential(double x_min, double x_max,
+                                               unsigned int Nx) {
   std::vector<double> values(Nx);
   std::vector<Point<dim>> eval_points(Nx);
 
   double Lx = x_max - x_min;
   double dx = Lx / Nx;
 
-  for(unsigned int i=0 ; i<Nx; ++i)
+  for (unsigned int i = 0; i < Nx; ++i)
     eval_points[i] = Point<1, double>(x_min + i * dx);
 
-  Utilities::MPI::RemotePointEvaluation<dim,dim> cache;
+  Utilities::MPI::RemotePointEvaluation<dim, dim> cache;
   cache.reinit(eval_points, triangulation, mapping);
 
   values = VectorTools::point_values<dim>(cache, dof_handler, solution);
@@ -175,45 +173,129 @@ std::vector<double> PoissonProblem<dim>::sample_electric_potential(
   return values;
 }
 
-// dealii Poisson
- 
-template<int dim>
-void PoissonProblem<dim>::create_mesh()
-{
+// // by GPT to re-re-re-check
+// template <int dim> std::vector<double> eval_solution_on_points(
+//         const std::vector<Vector<double>> &solutions,
+//         const unsigned int n,
+//         const std::vector<Point<dim>> &points, // need to be in [x_min, x_max]. I think....
+//         const std::vector<unsigned int> &cell_indices,
+//         const DoFHandler<dim> &dof_handler,
+//         const MappingQ<dim> &mapping)
+// {
+//     AssertIndexRange(n, solutions.size());
+//     Assert(points.size() == cell_indices.size(),
+//            ExcMessage("points and cell_indices must have same size"));
+//
+//     const Vector<double> &solution = solutions[n];
+//
+//     std::vector<double> result(points.size());
+//
+//     // Group points by cell (required for FEPointEvaluation efficiency)
+//     std::map<unsigned int, std::vector<unsigned int>> cell_to_point_ids;
+//
+//     for (unsigned int i = 0; i < points.size(); ++i)
+//         cell_to_point_ids[cell_indices[i]].push_back(i);
+//
+//     FEPointEvaluation<1, dim> evaluator(mapping,
+//                                         dof_handler.get_fe(),
+//                                         update_values);
+//
+//     std::vector<Point<dim>> cell_points;
+//     Vector<double> local_dofs(dof_handler.get_fe().dofs_per_cell);
+//
+//     for (const auto &entry : cell_to_point_ids)
+//     {
+//         const unsigned int cell_id = entry.first;
+//         const auto &point_ids = entry.second;
+//
+//         // these two lines bellow assume some order not sure how or why
+//         auto cell = dof_handler.begin_active();
+//         std::advance(cell, cell_id);
+//
+//         // extract points belonging to this cell
+//         cell_points.clear();
+//         cell_points.reserve(point_ids.size());
+//
+//         for (unsigned int id : point_ids)
+//             cell_points.push_back(points[id]);
+//
+//         std::vector<types::global_dof_index> indices(dof_handler.get_fe().n_dofs_per_cell());
+//         cell->get_dof_indices(indices);
+//
+//         for (unsigned int i=0;i<indices.size();++i)
+//           local_dofs[i] = solution[indices[i]];
+//
+//         // initialize evaluator on this cell
+//         evaluator.reinit(cell, cell_points);
+//
+//         evaluator.evaluate(local_dofs, EvaluationFlags::values);
+//
+//         for (unsigned int k = 0; k < point_ids.size(); ++k)
+//             result[point_ids[k]] = evaluator.get_value(k);
+//     }
+//
+//     return result;
+// }
 
-  GridGenerator::hyper_cube(triangulation,
-                            Parameters::X_DOMAIN_LEFT,
+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);
+}
+
+// dealii Poisson
+
+template <int dim> void PoissonProblem<dim>::create_mesh() {
+
+  GridGenerator::hyper_cube(triangulation, Parameters::X_DOMAIN_LEFT,
                             Parameters::X_DOMAIN_RIGHT);
 
+  std::vector<
+      GridTools::PeriodicFacePair<typename Triangulation<dim>::cell_iterator>>
+      periodic_faces;
 
-  std::vector<GridTools::PeriodicFacePair<
-      typename Triangulation<dim>::cell_iterator>> periodic_faces;
-
-  GridTools::collect_periodic_faces(triangulation,
-                                    0, 1,   // boundary IDs
-                                    0,      
-                                    periodic_faces);
+  GridTools::collect_periodic_faces(triangulation, 0, 1, // boundary IDs
+                                    0, periodic_faces);
 
   triangulation.add_periodicity(periodic_faces);
 
   triangulation.refine_global(Parameters::GLOBAL_REFINEMENT);
 }
 
-template <int dim>
-void PoissonProblem<dim>::setup_system()
-{
+template <int dim> void PoissonProblem<dim>::setup_system() {
 
   dof_handler.distribute_dofs(fe);
 
   constraints.clear();
-  
+
   DoFTools::make_hanging_node_constraints(dof_handler, constraints);
 
-  DoFTools::make_periodicity_constraints(dof_handler,
-                                         0, 1,
-                                         0,
-                                         constraints);
-  
+  DoFTools::make_periodicity_constraints(dof_handler, 0, 1, 0, constraints);
+
+  // Gauge fix for periodic Poisson:
+  // remove the constant nullspace by pinning one unconstrained DoF.
+  // (by Paul Wilhelm)
+  types::global_dof_index gauge_dof = numbers::invalid_dof_index;
+
+  for (types::global_dof_index i = 0; i < dof_handler.n_dofs(); ++i) {
+    if (!constraints.is_constrained(i)) {
+      gauge_dof = i;
+      break;
+    }
+
+
+  Assert(gauge_dof != numbers::invalid_dof_index,
+         ExcMessage("No unconstrained DoF found for gauge fixing."));
+
+  constraints.add_line(gauge_dof);
+  constraints.set_inhomogeneity(gauge_dof, 0.0);
+
   constraints.close();
 
   DynamicSparsityPattern dsp(dof_handler.n_dofs());
@@ -224,11 +306,21 @@ void PoissonProblem<dim>::setup_system()
 
   solution.reinit(dof_handler.n_dofs());
   system_rhs.reinit(dof_handler.n_dofs());
-}
 
+  fe_field_function =
+  std::make_unique<Functions::FEFieldFunction<dim>>(
+    dof_handler, solution, mapping);
+  }
+
+}
+/* (Mine)
 template <int dim>
 void PoissonProblem<dim>::assemble_system()
 {
+
+  system_matrix = 0;
+  system_rhs = 0;
+
   QGauss<dim>  quadrature_formula(fe.degree + 1);
   FEValues<dim> fe_values(fe, quadrature_formula,
                           update_values |
@@ -296,13 +388,56 @@ void PoissonProblem<dim>::assemble_system()
                                      solution,
                                      system_rhs);
 }
+ */
 
+// Paul's, mine's above
 
-template <int dim>
-void PoissonProblem<dim>::solve()
-{
+template <int dim> void PoissonProblem<dim>::assemble_system() {
+  system_matrix = 0;
+  system_rhs = 0;
 
-  SolverControl            solver_control(Parameters::CONVERGENCE_ITERATIONS, Parameters::CONVERGENCE_LIMIT);
+  QGauss<dim> quadrature_formula(fe.degree + 1);
+  FEValues<dim> fe_values(fe, quadrature_formula,
+                          update_values | update_gradients |
+                              update_quadrature_points | update_JxW_values);
+
+  const unsigned int dofs_per_cell = fe.n_dofs_per_cell();
+
+  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);
+
+  Assert(rhs_function != nullptr, ExcMessage("RHS function not set"));
+
+  for (const auto &cell : dof_handler.active_cell_iterators()) {
+    fe_values.reinit(cell);
+
+    cell_matrix = 0;
+    cell_rhs = 0;
+
+    for (const auto q : fe_values.quadrature_point_indices()) {
+      const double rho = rhs_function->value(fe_values.quadrature_point(q));
+
+      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->get_dof_indices(local_dof_indices);
+
+    constraints.distribute_local_to_global(
+        cell_matrix, cell_rhs, local_dof_indices, system_matrix, system_rhs);
+  }
+}
+
+template <int dim> void PoissonProblem<dim>::solve() {
+
+  SolverControl solver_control(Parameters::CONVERGENCE_ITERATIONS,
+                               Parameters::CONVERGENCE_LIMIT);
   SolverCG<Vector<double>> solver(solver_control);
 
   // PreconditionSSOR<SparseMatrix<double>> preconditioner;
@@ -310,28 +445,26 @@ void PoissonProblem<dim>::solve()
 
   // solver.solve(system_matrix, solution, system_rhs, preconditioner);
   solver.solve(system_matrix, solution, system_rhs, PreconditionIdentity());
-  // constraints.distribute(solution);
+  constraints.distribute(solution);
+
+
+  fe_field_function =
+  std::make_unique<Functions::FEFieldFunction<dim>>(
+    dof_handler, solution, mapping);
 }
 
-template <int dim>
-void PoissonProblem<dim>::initialize()
-{
-  create_mesh();      // build grid
-  setup_system();     // distribute DoFs and matrices
+template <int dim> void PoissonProblem<dim>::initialize() {
+  create_mesh();  // build grid
+  setup_system(); // distribute DoFs and matrices
 }
 
-template <int dim>
-void PoissonProblem<dim>::solve_step()
-{
+template <int dim> void PoissonProblem<dim>::solve_step() {
   assemble_system();
   solve();
 }
 
-
 // NuFI doesnt use this, kept only for testing PoissonProblem
-template <int dim>
-void PoissonProblem<dim>::run()
-{
+template <int dim> void PoissonProblem<dim>::run() {
   create_mesh();
   setup_system();
   assemble_system();
diff --git a/nufi/save_results.h b/nufi/save_results.h
index 2874edc..981158c 100644
--- a/nufi/save_results.h
+++ b/nufi/save_results.h
@@ -3,23 +3,24 @@
 
 #include <string>
 #include "nufi/nufi_solver.h"
+#include "nufi/poisson_problem.h"
 
 
 void save_f( const NuFISolver &solver,
     unsigned int n,
-                                    const double *E_coeffs,
+                                    const PoissonProblem<1> &poisson,
                                     unsigned int Nx_out,
                                     unsigned int Nv_out,
                                     const std::string &filename);
 
 void save_rho(const NuFISolver &solver,
     unsigned int n,
-                                    const double *E_coeffs,
+                                    const PoissonProblem<1> &poisson,
                                     unsigned int Nx_out,
                                     const std::string &filename);
 
 void save_Efield(unsigned int n,
-                                    const double *E_coeffs,
+                                    const PoissonProblem<1> &poisson,
                                     unsigned int Nx_out,
                                     const std::string &filename);
 
diff --git a/nufi/splines.h b/nufi/splines.h
deleted file mode 100644
index 73e8fd0..0000000
--- a/nufi/splines.h
+++ /dev/null
@@ -1,90 +0,0 @@
-#ifndef SPLINES_HP
-#define SPLINES_HP
-
-#include <cstddef>
-
-namespace splines1d
-{
-
-template <typename real>
-constexpr real faculty( size_t n ) noexcept
-{
-   return  (n > 1) ? real(n)*faculty<real>(n-1) : real(1);
-}
-
-template <typename real, size_t order, size_t derivative = 0>
-void N( real x, real *result, size_t stride = 1 ) noexcept
-{
-    static_assert( order > 0, "Splines must have order greater than zero." );
-    constexpr int n { order      };
-    constexpr int d { derivative };
-
-    if ( derivative >= order )
-        for ( size_t i = 0; i < order; ++i )
-            result[ i*stride ] = 0;
-
-    if ( n == 1 )
-    {
-        *result = 1;
-        return;
-    }
-
-    real v[n]; v[n-1] = 1;
-    for ( int k = 1; k < n - d; ++k )
-    {
-        v[n-k-1] = (1-x)*v[n-k];
-
-        for ( int i = 1-k; i < 0; ++i )
-            v[n-1+i] = (x-i)*v[n-1+i] + (k+1+i-x)*v[n+i];
-
-        v[n-1] *= x;
-    }
-
-    // Differentiate if necessary.
-    for ( size_t j = derivative; j-- > 0;  )
-    {
-        v[j] = -v[j+1];
-        for ( size_t i = j + 1; i < order - 1; ++i )
-            v[i] = v[i] - v[i+1];
-    }
-
-    constexpr real factor = real(1) / faculty<real>(order-derivative-1);
-    for ( size_t i = 0; i < order; ++i )
-        result[i*stride] = v[i]*factor;
-}
-
-template <typename real, size_t order, size_t derivative = 0>
-real eval( real x, const real *coefficients, size_t stride = 1 ) noexcept
-{
-    static_assert( order > 0, "Splines must have order greater than zero." );
-    static_assert( order > derivative, "Too high derivative requested." );
-    constexpr size_t n { order };
-    constexpr size_t d { derivative };
-
-    if ( d >= n ) return 0;
-    if ( n == 1 ) return *coefficients;
-
-    // Gather coefficients.
-    real c[ order ];
-    for ( size_t j = 0; j < order; ++j )
-        c[j] = coefficients[ stride * j ];
-
-    // Differentiate if necessary.
-    for ( size_t j = 1; j <= d; ++j )
-        for ( size_t i = n; i-- > j; )
-            c[i] = c[i] - c[i-1];
-
-    // Evaluate using de Boor’s algorithm.
-    for ( size_t j = 1; j < n-d; ++j )
-        for ( size_t i = n-d; i-- > j; )
-            c[d+i] = (x+n-d-1-i)*c[d+i] + (i-j+1-x)*c[d+i-1];
-
-    constexpr real factor = real(1) / faculty<real>(order-derivative-1);
-    return factor*c[n-1];
-}
-
-}
-
-#endif
-
-
diff --git a/rho_E.mp4 b/rho_E.mp4
index 2564d57..2730435 100644
Binary files a/rho_E.mp4 and b/rho_E.mp4 differ
diff --git a/src/blas.cc b/src/blas.cc
deleted file mode 100644
index 9107ecc..0000000
--- a/src/blas.cc
+++ /dev/null
@@ -1,95 +0,0 @@
-#include "nufi/blas.h"
-#include <cblas.h>
-
-namespace blas
-{
-
-double dot( const size_t n, const double *x, size_t incx,
-                            const double *y, size_t incy )
-{
-    return cblas_ddot(n,x,incx,y,incy);
-}
-
-float dot( const size_t n, const float *x, size_t incx,
-                           const float *y, size_t incy )
-{
-    return cblas_sdot(n,x,incx,y,incy);
-}
-
-void axpy( size_t n, double alpha, const double *x, size_t incx,
-                                         double *y, size_t incy )
-{
-    cblas_daxpy(n,alpha,x,incx,y,incy);
-}
-
-void axpy( size_t n, float  alpha, const float  *x, size_t incx,
-                                         float  *y, size_t incy )
-{
-    cblas_saxpy(n,alpha,x,incx,y,incy);
-}
-
-void scal( size_t n, double alpha, double *x, size_t incx )
-{
-    cblas_dscal(n,alpha,x,incx);
-}
-
-void scal( size_t n, float  alpha, float  *x, size_t incx )
-{
-    cblas_sscal(n,alpha,x,incx);
-}
-
-void copy( size_t n, const double *x, size_t incx, double *y, size_t incy )
-{
-    cblas_dcopy(n,x,incx,y,incy);
-}
-
-void copy( size_t n, const float  *x, size_t incx, float  *y, size_t incy )
-{
-    cblas_scopy(n,x,incx,y,incy);
-}
-
-void ger( const size_t M, const size_t N, const double alpha,
-          const double *X, const size_t incX, const double *Y, const size_t incY,
-          double *A, const size_t lda)
-{
-    cblas_dger( CblasColMajor, M, N, alpha, X, incX, Y, incY, A, lda );
-}
-
-void ger( const size_t M, const size_t N, const float alpha,
-          const float *X, const size_t incX, const float *Y, const size_t incY,
-          float *A, const size_t lda)
-{
-    cblas_sger( CblasColMajor, M, N, alpha, X, incX, Y, incY, A, lda );
-}
-
-void gemv( const char trans, size_t m, size_t n,
-		   double alpha, const double *a, size_t lda,
-           const double *x, size_t incx, double beta,
-           double *y, size_t incy )
-{
-    if ( trans == 'T' || trans == 'Y' )
-    {
-        cblas_dgemv( CblasColMajor, CblasTrans, m, n, alpha, a, lda, x, incx, beta, y, incy );
-    }
-    else
-    {
-        cblas_dgemv( CblasColMajor, CblasNoTrans, m, n, alpha, a, lda, x, incx, beta, y, incy );
-    }
-}
-
-void gemv( const char trans, size_t m, size_t n,
-		   float alpha, const float *a, size_t lda,
-           const float *x, size_t incx, float beta,
-           float *y, size_t incy )
-{
-    if ( trans == 'T' || trans == 'Y' )
-    {
-        cblas_sgemv( CblasColMajor, CblasTrans, m, n, alpha, a, lda, x, incx, beta, y, incy );
-    }
-    else
-    {
-        cblas_sgemv( CblasColMajor, CblasNoTrans, m, n, alpha, a, lda, x, incx, beta, y, incy );
-    }
-}
-
-}
diff --git a/src/nufi_solver.cc b/src/nufi_solver.cc
index 83ca6a1..deece5b 100644
--- a/src/nufi_solver.cc
+++ b/src/nufi_solver.cc
@@ -24,31 +24,24 @@ using namespace dealii;
 double NuFISolver::eval_ftilda(unsigned int n,
                                       double x,
                                       double u,
-                                      const double *E_coeffs) const
+                                      const PoissonProblem<1> &poisson) const
 {
   if ( n == 0 ) return f0(x,u);
   
-  const size_t order = Parameters::SPLINE_ORDER;
-  const size_t stride_x = 1;
-  const size_t stride_t = stride_x*(Nx + order - 1);
-
   double Ex;
-  const double *c;
 
   // We omit the initial half-step.
 
   while ( --n )
   {
       x  = x - Parameters::DT *u;
-      c  = E_coeffs + n*stride_t;
-      Ex = -eval<1>(x, c);
+      Ex = -eval(x, poisson);
       u  = u + Parameters::DT *Ex;
   }
 
   // The final half-step.
   x -= Parameters::DT*u;
-  c  = E_coeffs + n*stride_t;
-  Ex = -eval<1>(x, c);
+  Ex = -eval(x, poisson);
   u += 0.5*Parameters::DT*Ex;
 
   return f0(x,u);
@@ -57,34 +50,26 @@ double NuFISolver::eval_ftilda(unsigned int n,
 double NuFISolver::eval_f(unsigned int n,
                                       double x,
                                       double u,
-                                      const double *E_coeffs) const
+                                      const PoissonProblem<1> &poisson) const
 {
   if ( n == 0 ) return f0(x,u);
   
-  const size_t order = Parameters::SPLINE_ORDER;
-  const size_t stride_x = 1;
-  const size_t stride_t = stride_x*(Nx + order - 1);
-
   double Ex;
-  const double *c;
 
   // Initial half-step.
-  c  = E_coeffs + n*stride_t;
-  Ex = -eval<1>(x, c);
+  Ex = -eval(x, poisson);
   u += 0.5*Parameters::DT * Ex;
 
   while ( --n )
   {
       x  = x - Parameters::DT *u;
-      c  = E_coeffs + n*stride_t;
-      Ex = -eval<1>(x, c);
+      Ex = -eval(x, poisson);
       u  = u + Parameters::DT *Ex;
   }
 
   // The final half-step.
   x -= Parameters::DT*u;
-  c  = E_coeffs + n*stride_t;
-  Ex = -eval<1>(x, c);
+  Ex = -eval(x, poisson);
   u += 0.5*Parameters::DT*Ex;
 
   return f0(x,u);
@@ -92,7 +77,7 @@ double NuFISolver::eval_f(unsigned int n,
 
 double NuFISolver::eval_rho(const unsigned int n,
                           const double x,
-                          const double *E_coeffs,
+                          const PoissonProblem<1> &poisson,
                           const unsigned int Nv) const
 {
   const double dv = (Parameters::V_DOMAIN_RIGHT - Parameters::V_DOMAIN_LEFT) / Nv;
@@ -102,7 +87,7 @@ double NuFISolver::eval_rho(const unsigned int n,
 
 #pragma omp parallel for reduction (+ : integral)
   for (unsigned int i = 0; i < Nv; ++i)
-    integral += eval_ftilda(n, x, v_min + i * dv, E_coeffs);
+    integral += eval_ftilda(n, x, v_min + i * dv, poisson);
   
   return 1.0 - integral*dv;
 }
@@ -114,9 +99,6 @@ void NuFISolver::run()
   using std::abs;
   using std::max;
 
-  const size_t stride_t = Nx + order - 1;
-
-  std::unique_ptr<double[]> coeffs { new double[ Nt*stride_t ] {} };
   std::unique_ptr<double,decltype(std::free)*> rho { reinterpret_cast<double*>(std::aligned_alloc(64,sizeof(double)*Nx)), std::free };
 
   std::vector<double> int_E_squared;
@@ -138,13 +120,13 @@ void NuFISolver::run()
 
       // compute rho
 
-      double dx = Parameters::SPLINE_DX;
+      double dx = Parameters::CALC_DX;
       
       #pragma omp parallel for
     	for(size_t i = 0; i<Nx; i++)
     	{
         double x = Parameters::X_DOMAIN_LEFT + i*dx;
-        double ith_rho = eval_rho(it, x, coeffs.get(),Parameters::NV); 
+        double ith_rho = eval_rho(it, x, poisson, Parameters::NV); 
 
         AssertThrow(std::isfinite(ith_rho), ExcMessage("NaN detected in rho"));
     		rho.get()[i] =  ith_rho;
@@ -165,15 +147,12 @@ void NuFISolver::run()
 
 
       // interpolate and save current field
-      double* current_coeffs = coeffs.get() + it*stride_t;
-      interpolate<double, Parameters::SPLINE_ORDER>(current_coeffs, sampled_potential.data());
-
 
       std::vector<double> E_x(Nx,0.0) ;
       #pragma omp parallel for
       for(size_t ix=0; ix<Nx; ++ix)
       {
-        E_x[ix] = -eval<1>(Parameters::X_DOMAIN_LEFT+ix*dx, current_coeffs);
+        E_x[ix] = -eval(Parameters::X_DOMAIN_LEFT+ix*dx, poisson);
       }
 
       double timer_elapsed = timer.elapsed();
@@ -183,12 +162,12 @@ void NuFISolver::run()
       if (it % Parameters::PLOT_FREQUENCY == 0)
       {
         std::cout << "Saving results...   ";
-        save_f(*this, it, coeffs.get(), Parameters::SPLINE_NX, Parameters::NV, "results/ftilda_" + std::to_string(it) + ".dat");
-        save_rho(*this, it, coeffs.get(), Parameters::SPLINE_NX, "results/rho_" + std::to_string(it) + ".dat");
+        save_f(*this, it, poisson, Parameters::PLOT_NX, Parameters::NV, "results/ftilda_" + std::to_string(it) + ".dat");
+        save_rho(*this, it, poisson, Parameters::PLOT_NX, "results/rho_" + std::to_string(it) + ".dat");
         // save_Efield(it, coeffs.get(), 128, "results/field_" + std::to_string(it) + ".dat");
         save_space_vector(E_x, "field", it);
 
-        double int_val = 0.5 * integral_space_vector_squared(current_coeffs);
+        double int_val = 0.5 * integral_space_vector_squared(poisson);
         int_E_squared.push_back(int_val);  
         save_space_vector(int_E_squared, "electricint", it);
         std::cout << "Time since start = "<< total_time<<"\n\n";
diff --git a/src/save_results.cc b/src/save_results.cc
index 8c6d063..6d1b31d 100644
--- a/src/save_results.cc
+++ b/src/save_results.cc
@@ -6,11 +6,13 @@
 #include <string>
 #include <vector>
 #include "nufi/nufi_solver.h"
+#include "nufi/poisson_problem.h"
+#include "nufi/fields.h"
 
 
 void save_f( const NuFISolver &solver,
     unsigned int n,
-                                    const double *E_coeffs,
+                                    const PoissonProblem<1> &poisson,
                                     unsigned int Nx_out,
                                     unsigned int Nv_out,
                                     const std::string &filename)
@@ -38,7 +40,7 @@ void save_f( const NuFISolver &solver,
       {
           double v = vmin + (j + 0.5)*dv;
 
-          double val = solver.eval_f(n, x, v, E_coeffs);
+          double val = solver.eval_f(n, x, v, poisson);
 
           file << val;
 
@@ -54,7 +56,7 @@ void save_f( const NuFISolver &solver,
 
 void save_rho(const NuFISolver &solver,
     unsigned int n,
-                                    const double *E_coeffs,
+                                    const PoissonProblem<1> &poisson,
                                     unsigned int Nx_out,
                                     const std::string &filename)
 {
@@ -69,15 +71,15 @@ void save_rho(const NuFISolver &solver,
 
   for (unsigned int i = 0; i < Nx_out; ++i, xmin += dx)
   {
-      double val = solver.eval_rho(n, xmin, E_coeffs);
+      double val = solver.eval_rho(n, xmin, poisson);
       file << val;
       file << "\n";
   }
   file.close();
 }
 
-void save_Efield(unsigned int n,
-                                    const double *E_coeffs,
+void save_Efield([[maybe_unused]]unsigned int n,
+                                    const PoissonProblem<1> &poisson,
                                     unsigned int Nx_out,
                                     const std::string &filename)
 {
@@ -88,17 +90,13 @@ void save_Efield(unsigned int n,
   double dx = (xmax - xmin) / Nx_out;
 
   // select from E_coeffs
-  const size_t stride_x = 1;
-  const size_t stride_t = stride_x*(Parameters::SPLINE_NX + Parameters::SPLINE_ORDER - 1);
-  const double *c;
-  c  = E_coeffs + n*stride_t;
 
   file << Nx_out << "\n";
   file << xmin << " " << xmax << "\n";
 
   for (unsigned int i = 0; i < Nx_out; ++i, xmin += dx)
   {
-      double val = -eval<1>(xmin, c);
+      double val = -eval(xmin, poisson);
       file << val;
       file << "\n";
   }