diff --git a/CMakeLists.txt b/CMakeLists.txt
index 2732c3c..60ddc93 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -21,6 +21,7 @@ deal_ii_initialize_cached_variables()
 add_library(nufi_lib
   src/nufi_solver.cc
   src/save_results.cc
+  src/blas.cc
 )
 
 target_include_directories(nufi_lib PUBLIC
diff --git a/Makefile b/Makefile
index 38f5543..1b6e392 100644
--- a/Makefile
+++ b/Makefile
@@ -142,6 +142,30 @@ 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
 
@@ -224,6 +248,9 @@ 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 ccc2bf6..2c21e05 100644
Binary files a/libnufi_lib.a and b/libnufi_lib.a differ
diff --git a/nufi/fields.h b/nufi/fields.h
index 74cf3cc..86f782d 100644
--- a/nufi/fields.h
+++ b/nufi/fields.h
@@ -65,14 +65,14 @@ double eval(double x, const double *coeffs) noexcept
 }
 
 template <typename real, size_t order>
-void interpolate( real *coeffs, const real *values) // Least Squares needs to be made
+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 // STRUCT AND CONFIG NEEDS TO BE REVIEWED 
+    struct mat_t
     {
         real  N[ order ];
 
@@ -101,7 +101,7 @@ void interpolate( real *coeffs, const real *values) // Least Squares needs to be
         }
     };
 
-    struct transposed_mat_t // STRUCT AND CONFIG NEEDS TO BE REVIEWED 
+    struct transposed_mat_t
     {
         real  N[ order ];
 
diff --git a/nufi/nufi_solver.h b/nufi/nufi_solver.h
index 2a69624..756cf2b 100644
--- a/nufi/nufi_solver.h
+++ b/nufi/nufi_solver.h
@@ -9,7 +9,7 @@
 
 #include "nufi/parameters.h"
 #include "nufi/poisson_problem.h"
-#include "nufi/fields.h"
+#include "nufi/fields.h" //dont remove
 
 using namespace dealii;
 
diff --git a/nufi/parameters.h b/nufi/parameters.h
index ebb11a3..93f9b02 100644
--- a/nufi/parameters.h
+++ b/nufi/parameters.h
@@ -15,10 +15,10 @@ namespace Parameters
   constexpr double V_DOMAIN_LEFT = -10.0;
   constexpr double V_DOMAIN_RIGHT = 10.0;
 
-  constexpr unsigned int NV = 512;
+  constexpr unsigned int NV = 1024;
 
   constexpr unsigned int GLOBAL_REFINEMENT = 8;
-  constexpr unsigned int FE_DEGREE = 4;
+  constexpr unsigned int FE_DEGREE = 3;
 
   constexpr double EPS = 0.01;
   constexpr double WAVE_NR = 0.5;
@@ -30,7 +30,7 @@ namespace Parameters
 
 
   //spline options
-  constexpr int SPLINE_NX = 512;
+  constexpr int SPLINE_NX = 1024;
   constexpr double SPLINE_DX = LX/SPLINE_NX;
   constexpr size_t SPLINE_ORDER = 4;
 
diff --git a/nufi/poisson_problem.h b/nufi/poisson_problem.h
index 0f5bdfc..f118050 100644
--- a/nufi/poisson_problem.h
+++ b/nufi/poisson_problem.h
@@ -30,6 +30,7 @@
 
 #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 <memory>
@@ -58,8 +59,7 @@ 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(const PoissonProblem<dim> &problem,  // sampling to save as spline
-                                            unsigned int Nx,
+  std::vector<double> sample_electric_field(unsigned int Nx,
                                             double x_min,
                                             double x_max);
 
@@ -103,36 +103,37 @@ PoissonProblem<dim>::PoissonProblem(unsigned int degree)
 
 template <int dim>
 std::vector<double> PoissonProblem<dim>::sample_electric_field(
-    const PoissonProblem<dim> &problem,
     unsigned int Nx,
     double x_min,
     double x_max)
 {
+  // const auto &dof_handler = problem.get_dof_handler();
+  // const auto &solution = problem.get_solution();
 
-    const auto &dof_handler = problem.get_dof_handler();
-    const auto &solution = problem.get_solution();
+  this -> get_solution();
+  this -> get_dof_handler();
 
-    Functions::FEFieldFunction<dim, Vector<double>>
-      field_function(dof_handler, solution, mapping);
+  Functions::FEFieldFunction<dim, Vector<double>>
+    field_function(dof_handler, solution, mapping);
 
-    std::vector<double> values(Nx);
+  std::vector<double> values(Nx);
 
-    double Lx = x_max - x_min;
-    double dx = Lx / Nx;
+  double Lx = x_max - x_min;
+  double dx = Lx / Nx;
 
-    for (unsigned int i = 0; i < Nx; ++i)
-    {
-        double x = x_min + i * dx;
+  for (unsigned int i = 0; i < Nx; ++i)
+  {
+      double x = x_min + i * dx;
 
-        Point<dim> p;
-        p[0] = x;
+      Point<dim> p;
+      p[0] = x;
 
-        Tensor<1, dim> grad = field_function.gradient(p);
+      Tensor<1, dim> grad = field_function.gradient(p);
 
-        values[i] = -grad[0];  // E = -dφ/dx
-    }
+      values[i] = -grad[0];  // E = -dφ/dx
+  }
 
-    return values;
+  return values;
 }
 
 // dealii Poisson
@@ -145,19 +146,23 @@ void PoissonProblem<dim>::create_mesh()
                             Parameters::X_DOMAIN_LEFT,
                             Parameters::X_DOMAIN_RIGHT);
 
+  // TO CHECK:
+  //
   // Make x-dim boundaries periodic
-  Tensor<1, dim> offset;
-  std::vector<GridTools::PeriodicFacePair<
-              typename Triangulation<dim>::cell_iterator>> periodicity_vector;
-
-  GridTools::collect_periodic_faces(triangulation,
-                                    0,
-                                    1,
-                                    0,
-                                    periodicity_vector,
-                                    offset);
-  
-  triangulation.add_periodicity(periodicity_vector);
+  // Tensor<1, dim> offset;
+  // std::vector<GridTools::PeriodicFacePair<
+  //             typename Triangulation<dim>::cell_iterator>> periodicity_vector;
+  //
+  // GridTools::collect_periodic_faces(triangulation,
+  //                                   0,
+  //                                   1,
+  //                                   0,
+  //                                   periodicity_vector,
+  //                                   offset);
+  //
+  // triangulation.add_periodicity(periodicity_vector);
+  //
+  // END CHECK
 
   triangulation.refine_global(Parameters::GLOBAL_REFINEMENT);
 }
@@ -168,50 +173,29 @@ void PoissonProblem<dim>::setup_system()
 
   dof_handler.distribute_dofs(fe);
 
-  constraints.clear();
-  DoFTools::make_hanging_node_constraints(dof_handler, constraints);
-
-  // 'boundary' condition phi(x_0) = 0 
-  constraints.add_line(0);
-  constraints.set_inhomogeneity(0, 0.0);
-
-  constraints.close();
+  // TO CHECK:
+  //
+  // constraints.clear();
+  // DoFTools::make_hanging_node_constraints(dof_handler, constraints);
+  //
+  // // 'boundary' condition phi(x_0) = 0 
+  // constraints.add_line(0);
+  // constraints.set_inhomogeneity(0, 0.0);
+  //
+  // constraints.close();
+  //
+  // END CHECK
 
   DynamicSparsityPattern dsp(dof_handler.n_dofs());
   DoFTools::make_sparsity_pattern(dof_handler, dsp, constraints);
   sparsity_pattern.copy_from(dsp);
 
   system_matrix.reinit(sparsity_pattern);
+
   solution.reinit(dof_handler.n_dofs());
   system_rhs.reinit(dof_handler.n_dofs());
 }
 
-// =-=-=-=-= E_field = -dPhi/dx =-=-=-=-=
-
-template <int dim>
-class ElectricFieldPostprocessor : public DataPostprocessorVector<dim>
-{
-public:
-  ElectricFieldPostprocessor()
-    : DataPostprocessorVector<dim>("electric_field", update_gradients)
-  {}
-
-  virtual void evaluate_scalar_field(
-    const DataPostprocessorInputs::Scalar<dim> &input_data,
-    std::vector<Vector<double>> &computed_quantities) const override
-  {
-    AssertDimension(input_data.solution_gradients.size(),
-                    computed_quantities.size());
-
-    for (unsigned int p = 0; p < input_data.solution_gradients.size(); ++p)
-      {
-        AssertDimension(computed_quantities[p].size(), dim);
-        for (unsigned int d = 0; d < dim; ++d)
-          computed_quantities[p][d] = -input_data.solution_gradients[p][d];
-      }
-  }
-};
-
 template <int dim>
 void PoissonProblem<dim>::assemble_system()
 {
@@ -223,7 +207,6 @@ void PoissonProblem<dim>::assemble_system()
                           update_JxW_values);
 
   const unsigned int dofs_per_cell = fe.n_dofs_per_cell();
-  const unsigned int n_q_points    = quadrature_formula.size();
 
   FullMatrix<double> cell_matrix(dofs_per_cell, dofs_per_cell);
   Vector<double>     cell_rhs(dofs_per_cell);
@@ -234,35 +217,54 @@ void PoissonProblem<dim>::assemble_system()
   for (const auto &cell : dof_handler.active_cell_iterators())
   {
     fe_values.reinit(cell);
+
     cell_matrix = 0;
     cell_rhs    = 0;
 
-    for (unsigned int q = 0; q < n_q_points; ++q)
+    for (const auto q : fe_values.quadrature_point_indices())
     {
       const double rho = rhs_function->value(fe_values.quadrature_point(q));
 
-      for (unsigned int i = 0; i < dofs_per_cell; ++i)
-      {
-        for (unsigned int j = 0; j < dofs_per_cell; ++j)
+      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);
+            (fe_values.shape_grad(i, q) * // grad phi_i(x_q)
+             fe_values.shape_grad(j, q) * // grad phi_j(x_q)
+             fe_values.JxW(q));           // dx
+
+      for (const unsigned int i : fe_values.dof_indices())
+        cell_rhs(i) += (fe_values.shape_value(i, q) * // phi_i(x_q)
+                        rho *                         // f(x_q)
+                        fe_values.JxW(q));            // dx
+
 
-        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);
+    // constraints.distribute_local_to_global(cell_matrix,
+    //                                        cell_rhs,
+    //                                        local_dof_indices,
+    //                                        system_matrix,
+    //                                        system_rhs);
+    for (const unsigned int i : fe_values.dof_indices())
+      for (const unsigned int j : fe_values.dof_indices())
+        system_matrix.add(local_dof_indices[i],
+                          local_dof_indices[j],
+                          cell_matrix(i, j));
+
+    for (const unsigned int i : fe_values.dof_indices())
+      system_rhs(local_dof_indices[i]) += cell_rhs(i);
+
   }
+  std::map<types::global_dof_index, double> boundary_values;
+  VectorTools::interpolate_boundary_values(dof_handler,
+                                           types::boundary_id(0),
+                                           Functions::ZeroFunction<1>(),
+                                           boundary_values);
+  MatrixTools::apply_boundary_values(boundary_values,
+                                     system_matrix,
+                                     solution,
+                                     system_rhs);
 }
 
 
@@ -270,14 +272,15 @@ template <int dim>
 void PoissonProblem<dim>::solve()
 {
 
-  SolverControl            solver_control(1000, 1e-12);
+  SolverControl            solver_control(5000, 1e-12);
   SolverCG<Vector<double>> solver(solver_control);
 
-  PreconditionSSOR<SparseMatrix<double>> preconditioner;
-  preconditioner.initialize(system_matrix, 1.2);
+  // PreconditionSSOR<SparseMatrix<double>> preconditioner;
+  // preconditioner.initialize(system_matrix, 1.2);
 
-  solver.solve(system_matrix, solution, system_rhs, preconditioner);
-  constraints.distribute(solution);
+  // solver.solve(system_matrix, solution, system_rhs, preconditioner);
+  solver.solve(system_matrix, solution, system_rhs, PreconditionIdentity());
+  // constraints.distribute(solution);
 }
 
 template <int dim>
diff --git a/src/blas.cc b/src/blas.cc
new file mode 100644
index 0000000..9107ecc
--- /dev/null
+++ b/src/blas.cc
@@ -0,0 +1,95 @@
+#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 4be1506..2bb11e0 100644
--- a/src/nufi_solver.cc
+++ b/src/nufi_solver.cc
@@ -1,5 +1,6 @@
 #include "nufi/nufi_solver.h"
 
+#include <algorithm>
 #include <cmath>
 #include <cstdlib>
 #include <deal.II/base/point.h>
@@ -10,6 +11,7 @@
 #include <memory>
 #include <ostream>
 #include <cstddef>
+#include <vector>
 
 #include "nufi/parameters.h"
 #include "nufi/save_results.h"
@@ -90,32 +92,19 @@ void NuFISolver::run()
 
     	for(size_t i = 0; i<Nx; i++, x+=dx)
     	{
-        double ith_rho = eval_rho(it, x, coeffs.get(), Parameters::NV); 
+        double ith_rho = eval_rho(it, x, coeffs.get(),Parameters::NV); 
         AssertThrow(std::isfinite(ith_rho), ExcMessage("NaN detected in rho"));
     		rho.get()[i] =  ith_rho;
     	}
 
-      // bool rho_written = false;
-      //
-      // if (!rho_written)
-      // {
-      //     std::ofstream rho_file("rho_initial.dat");
-      //
-      //     if (!rho_file)
-      //         throw std::runtime_error("Could not open rho_initial.dat");
-      //
-      //     rho_file.precision(16);
-      //     rho_file << std::scientific;
-      //
-      //     for (size_t i = 0; i < Nx; ++i)
-      //         rho_file << rho.get()[i] << "\n";
-      //
-      //     rho_written = true;
-      // }
-
       poisson.set_rhs_function(std::make_unique<ChargeDensity_NuFI<1>>(rho.get(), Parameters::SPLINE_NX));
       poisson.solve_step();
 
+      std::vector<double> E_vals = poisson.sample_electric_field(Parameters::SPLINE_NX, Parameters::X_DOMAIN_LEFT, Parameters::X_DOMAIN_RIGHT);
+
+      double* current_coeffs = coeffs.get() + it*stride_t;
+      interpolate<double, Parameters::SPLINE_ORDER>(current_coeffs, E_vals.data());
+
       if (it % Parameters::PLOT_FREQUENCY == 0)
       {
         std::cout << "Saving results... \n\n";