diff --git a/poisson2d/README.md b/poisson2d/README.md
index 14ae309..086c4b4 100644
--- a/poisson2d/README.md
+++ b/poisson2d/README.md
@@ -15,17 +15,21 @@ The grid is a 2-dimensional array of size MxM. The timings for different values
 
 | Language            | M=100           | M=200                 | M=300                  |
 |---------------------|-----------------|-----------------------|------------------------|
-| Python (pure)       | 276             | n/a                   | n/a                    |
-| Cython              | 1.02            | 32.8                  | 229                    |
-| Fortran (naive)     | 0.34            | 13.2                  | 69.7                   |
-| Fortran (optimized) | 0.18            | 6.25                  | 31.4                   |
-| C (naive)           | 0.42*           | 7.25                  | 33.7                   |
-| C (optimized)       | 0.37*           | 6.80                  | 32.8                   |
-
-* For all of these results, the amount of iterations performed by the respective codes was approximately
-the same, with the exception of the 100x100 grid C codes, which did nearly a double amount of iterations
-compared to the rest, for some reason. The timings are on AMD Ryzen 5 3600 @3.6GHz, using WSL2 on Windows 10.
+| Python (pure)       | 198             | n/a                   | n/a                    |
+| Cython              |	0.46            | 5.77                  | 43.9                   |
+| Fortran (trad1)     | 0.23            | 3.59                  | 18.3                   |
+| Fortran (trad2)     | 0.19            | 4.29                  | 19.6                   |
+| C (naive)           | 0.18            | 2.46                  | 11.6                   |
+| C (optimized)       | 0.13            | 1.83                  | 8.26                   |
+| Fortran (modern)    | 0.19            | n/a                   | 8.08                   |
+| Python (vectorized) | 0.92            | 10.8                  | 65.7                   |
+
+* For all of these results, the amount of iterations performed was exactly the same.
+ The timings are on AMD Ryzen 5 3600 @3.6GHz, using WSL2 on Windows 10.
 
 Some thoughts on the code at https://runningcrocodile.fi/articles/pythonperformance.html . Also, there was
 discussion about this problem at Discourse: https://fortran-lang.discourse.group/t/performance-c-vs-fortran/1461 
 with important contributions from numerous users (Mason, Beliavsky, septc, implicitall, nncarlson, han190 and pmk)
+
+Codes written by a running crocodile (traditional1.f90, traditional2.f90, naive.c, optimized.c, poisson.py, poissonvectorized.py, poisson.pyx)
+and Damian Rouson (modern.f90) with the aforementioned help from Discourse.
diff --git a/poisson2d/modern.f90 b/poisson2d/modern.f90
new file mode 100644
index 0000000..81e00fd
--- /dev/null
+++ b/poisson2d/modern.f90
@@ -0,0 +1,97 @@
+module modern_func
+  !! Code by Damian Rouson, @rouson.
+  !! Modified by ARC, @arunningcroc
+
+  !! Define a scalar function over 2-space.
+  implicit none
+
+  private
+  public :: dp, run_modern
+
+  integer, parameter :: precision = 15, range = 50
+  integer, parameter :: dp = selected_real_kind(precision, range)
+contains
+  pure real(dp) function rho(x,y)
+    real(dp), intent(in) :: x,y
+    rho = -(6.0_dp*x*y*(1.0_dp-y)-2.0_dp*x**3.0_dp)
+  end function
+
+  pure elemental real(dp) function boundary(y)
+    real(dp),intent(in)   :: y
+    boundary = y*(1-y)
+  end function
+
+  pure real(dp) function analytical_sol(x,y)
+    real(dp), intent(in) :: x,y
+    analytical_sol = y*(1-y)*x**3
+  end function
+
+  integer function run_modern(M)
+    real(dp)    :: dx, rho_sampled(M,M)
+    integer     :: M,i,j
+
+    dx = 1.0/(M-1)
+    associate( dy => (dx) ) ! Associating with an expression provides immutable state so dy cannot be inadvertently redefined.
+
+      do concurrent(i=1:M, j=1:M)
+        rho_sampled(i,j) = rho(i*dx,j*dy)
+      end do
+   
+      block ! Tighten the scoping to declutter the code above.
+        real(dp)                         :: delta_phi, avgerror
+        real(dp), parameter              ::  tolerance=1E-8_dp, analytical_tol = 5.0e-3_dp
+        real(dp), dimension(0:M-1,0:M-1) :: phi_prime, phi, analytical
+        integer                          :: iteration
+  
+        phi = 0.
+        phi(M-1,:) = boundary([(i*dx,i=0,M-1)])
+
+        do concurrent(i=0:M-1, j=0:M-1)
+          analytical(i,j) = analytical_sol(i*dx,j*dx)
+        end do
+
+        phi_prime([0,M-1], 1:M-2) = phi([0,M-1], 1:M-2) ! Initialize only boundary values except corners. (Internal values will
+        phi_prime(1:M-2, [0,M-1]) = phi(1:M-2, [0,M-1]) ! be overwritten in the first iteration. Corners will never be used.)
+        
+        delta_phi = 1.0_dp !tolerance + epsilon(tolerance) ! Ensure at least 1 iteration.
+        iteration = 0
+        do while (delta_phi > tolerance )
+          iteration = iteration + 1
+          do concurrent(i=1:M-2, j=1:M-2) ! Compute updated solution estimate at internal points.
+            phi_prime(i,j) =   (phi(i+1,j) + phi(i-1,j) + phi(i,j+1) + phi(i,j-1))/4._dp &
+                             + (dx/2._dp)*(dy/2._dp)*rho_sampled(i,j)
+          end do
+          delta_phi = maxval(abs(phi_prime - phi))
+          phi(1:M-2, 1:M-2) = phi_prime(1:M-2, 1:M-2) ! Update internal values.
+          !print *, iteration
+        end do
+        phi = phi/maxval(abs(phi))
+        analytical = analytical/maxval(abs(analytical))
+        avgerror = sum(abs(phi-analytical))/(M*M)
+        if(avgerror > analytical_tol) then
+            print *, "Failed to match analytical solution", avgerror
+        end if
+        run_modern = iteration
+      end block
+      
+
+    end associate
+  end function
+end module 
+
+
+program poisson
+  !! Solve the 2D Poisson equation using a smoothing operation to iterate.
+  use modern_func, only: run_modern, dp
+  implicit none
+
+  integer             :: iter
+  integer, parameter  :: M=170
+  real(dp)            :: t_start, t_end
+
+  
+  call cpu_time(t_start)
+  iter = run_modern(M)
+  call cpu_time(t_end)
+  print *, t_end-t_start, iter
+end program
\ No newline at end of file
diff --git a/poisson2d/naive.c b/poisson2d/naive.c
index 9d72f22..6ed2091 100644
--- a/poisson2d/naive.c
+++ b/poisson2d/naive.c
@@ -3,7 +3,14 @@
 #include <string.h>
 #include <math.h>
 #include <time.h>
-#define M 300
+#define M 100
+#define analytical_tol 5.0e-3
+
+/*
+ *
+ * Code by ARC, @arunningcroc
+ *
+ */
 void swap(double (*a)[M], double (*b)[M], int n)
 {
   double swp;
@@ -28,38 +35,69 @@ double mmax(double (*phi)[M], double (*phip)[M], int n)
   }
   return max;
 }
+void normalize(double (*phi)[M]) {
+  double max = -50000;
+  for(int i=0; i<M; i++) {
+    for(int j=0; j<M; j++) {
+      if(fabs(phi[i][j]) > max) max = fabs(phi[i][j]);
+    }
+  }
+  for(int i=0; i<M; i++) {
+    for(int j=0; j<M; j++) {
+      phi[i][j] = phi[i][j]/max;
+    }
+  }
+  
+}
 double rho(double x, double y)
 {
-  double s1 = 0.6;
-  double e1 = 0.8;
-	double s2 = 0.2;
-	double e2 = 0.4;
-
-	if (x > s1 && x < e1 && y > s1 && y < e1) {
-		return 1.0;
-	} else if (x > s2 && x < e2 && y > s2 && y < e2 ) {
-		return -1.0;
-	} else {
-		return 0.0;
-	}
+  return -(6.0*x*y*(1.0-y)-2.0*pow(x,3.0));
 }
-void run(double toler, double a)
+double boundary(double y)
+{
+  return y*(1.0-y);
+}
+double analytical_sol(double x,double y)
+{
+  return y*(1.0-y)*pow(x,3.0);
+}
+double get_average_error(double (*sol)[M], double (*analytical)[M]) {
+  double sum = 0.0;
+  normalize(sol);
+  normalize(analytical);
+  for(int i=0; i<M; i++) {
+    for(int j=0; j<M; j++) {
+      sum += fabs(sol[i][j]-analytical[i][j]);
+    }
+  }
+  return sum/(M*M);
+}
+int run(double toler, double a)
 {
-  double epsilon0 = 8.85e-12;
   double a2;
 
   double (*phi)[M];
   double (*phip)[M]; 
-  double (*rhoa)[M]; 
+  double (*rhoa)[M];
+  double (*analytical)[M];
   phi = malloc(sizeof(double[M][M]));
   phip = malloc(sizeof(double[M][M]));
   rhoa = malloc(sizeof(double[M][M]));
-
+  analytical = malloc(sizeof(double[M][M]));
+  for(int i=0; i<M; i++) {
+    for (int j=0; j<M; j++) {
+      analytical[i][j] = analytical_sol(i*a,j*a);
+    }
+  }
   int iter = 0;
 
   memset(phip, 0, sizeof(phip[0][0])*M*M);
   memset(phi, 0, sizeof(phi[0][0])*M*M);
   
+  for (int j=0; j < M; j++) {
+    phip[M-1][j] = j*a*(1-j*a);
+    phi[M-1][j] = j*a*(1-j*a);
+  }
   double delta = 1.0;
   a2 = pow(a,2.0);
   while (delta > toler) {
@@ -69,22 +107,25 @@ void run(double toler, double a)
       for (int j=1; j < M-1; j++) {
         phip[i][j] = (phi[i+1][j] + phi[i-1][j] +
                       phi[i][j+1] + phi[i][j-1])/4.0 +
-                      a2/(4.0 * epsilon0)*rho(i*a,j*a);
+                      a2*rho(i*a,j*a)/4.0;
       }
     }
     delta = mmax(phi, phip, M);
     swap(phi, phip, M);
 
   }
-  printf("iters %d", iter);
-
+  if(get_average_error(phi, analytical) > analytical_tol) {
+    printf("Failed to reach analytical solution, error %f \n", get_average_error(phi, analytical));
+  }
+  return iter;
 }
 
 int main(int argc, char *argv[])
 {
+  int iter;
   clock_t start = clock();
-  run(1e-6, 0.01);
+  iter = run(1e-8, 1.0/(M-1));
   clock_t end = clock();
   double total = ((double)(end - start)) / CLOCKS_PER_SEC;
-  printf("Execution time: %f\n",total);
+  printf("Execution time %f, iters: %d\n",total,iter);
 }
diff --git a/poisson2d/naive.f90 b/poisson2d/naive.f90
deleted file mode 100644
index 3c34692..0000000
--- a/poisson2d/naive.f90
+++ /dev/null
@@ -1,51 +0,0 @@
-module rhofunc
-implicit none
-public
-integer, parameter :: dp=kind(0.d0)
-contains
-    pure real(dp) function rho(x,y)
-        real(dp), intent(in) :: x,y
-        if (x > 0.6_dp .and. x < 0.8_dp .and. y > 0.6_dp .and. y<0.8) then
-            rho = 1.0_dp
-        else if (x> 0.2_dp .and. x<0.4_dp .and. y>0.2_dp .and. y<0.4_dp) then
-            rho = -1.0_dp
-        else
-            rho = 0.0_dp
-        end if
-    end function
-
-end module
-
-program poisson
-use rhofunc, only: rho
-implicit none
-integer, parameter :: dp=kind(0.d0), M=300
-integer            :: i,j, iter
-real(dp),parameter :: epsilon0=8.85E-12_dp, target=1E-6_dp, a=0.01
-real(dp)           :: delta, b, e, phiprime(M,M), phi(M,M), a2, rhoarr(M,M), temp(M,M)
-
-
-delta = 1.0_dp
-iter = 0
-call cpu_time(b)
-phiprime(:,:) = 0.0_dp
-phi(:,:) = 0.0_dp
-
-do while (delta > target )
-    iter = iter + 1
-    a2 = a**2.0_dp
-    do i=2, M-1
-        do j=2, M-1
-            phiprime(i,j) = (phi(i+1,j) + phi(i-1,j) + phi(i,j+1) + phi(i,j-1))/4.0_dp &
-            + a2/4.0_dp/epsilon0*rho(i*a,j*a)
-        end do
-    end do
-    delta = maxval(abs(phiprime - phi))
-    temp = phi
-    phi = phiprime 
-    phiprime = temp
-
-end do
-call cpu_time(e)
-print *, e-b, iter
-end program
\ No newline at end of file
diff --git a/poisson2d/optimized.c b/poisson2d/optimized.c
index 0dddfec..b758a1f 100644
--- a/poisson2d/optimized.c
+++ b/poisson2d/optimized.c
@@ -3,7 +3,14 @@
 #include <string.h>
 #include <math.h>
 #include <time.h>
-#define M 300
+#define M 200
+#define analytical_tol 5.0e-3
+
+/*
+ *
+ * Code by ARC, @arunningcroc
+ *
+ */
 void swap(double (*a)[M], double (*b)[M], int n)
 {
   double swp;
@@ -26,33 +33,60 @@ double mmax(double (*phi)[M], double (*phip)[M], int n)
   }
   return max;
 }
+void normalize(double (*phi)[M]) {
+  double max = -50000;
+  for(int i=0; i<M; i++) {
+    for(int j=0; j<M; j++) {
+      if(phi[i][j] > max) max = phi[i][j];
+    }
+  }
+  for(int i=0; i<M; i++) {
+    for(int j=0; j<M; j++) {
+      phi[i][j] /= max;
+    }
+  }
+  
+}
+double get_average_error(double (*sol)[M], double (*analytical)[M]) {
+  double sum = 0.0;
+  normalize(sol);
+  normalize(analytical);
+  for(int i=0; i<M; i++) {
+    for(int j=0; j<M; j++) {
+      sum += fabs(sol[i][j]-analytical[i][j]);
+    }
+  }
+  return sum/(M*M);
+}
 double rho(double x, double y)
 {
-  double s1 = 0.6;
-  double e1 = 0.8;
-	double s2 = 0.2;
-	double e2 = 0.4;
-
-	if (x > s1 && x < e1 && y > s1 && y < e1) {
-		return 1.0;
-	} else if (x > s2 && x < e2 && y > s2 && y < e2 ) {
-		return -1.0;
-	} else {
-		return 0.0;
-	}
+  return -(6.0*x*y*(1.0-y)-2.0*pow(x,3.0));
 }
-void run(double toler, double a)
+double boundary(double y)
+{
+  return y*(1.0-y);
+}
+double analytical_sol(double x,double y)
+{
+  return y*(1.0-y)*pow(x,3.0);
+}
+int run(double toler, double a)
 {
-  double epsilon0 = 8.85e-12;
   double a2;
 
   double (*phi)[M];
   double (*phip)[M]; 
   double (*rhoa)[M]; 
+  double (*analytical)[M];
   phi = malloc(sizeof(double[M][M]));
   phip = malloc(sizeof(double[M][M]));
   rhoa = malloc(sizeof(double[M][M]));
-
+  analytical = malloc(sizeof(double[M][M]));
+  for(int i=0; i<M; i++) {
+    for (int j=0; j<M; j++) {
+      analytical[i][j] = analytical_sol(i*a,j*a);
+    }
+  }
   int iter = 0;
 
   memset(phip, 0, sizeof(phip[0][0])*M*M);
@@ -62,6 +96,10 @@ void run(double toler, double a)
       rhoa[i][j] =  rho(i*a,j*a);
     }
   }
+  for (int j=0; j < M; j++) {
+    phip[M-1][j] = j*a*(1-j*a);
+    phi[M-1][j] = j*a*(1-j*a);
+  }
   double delta = 1.0;
   a2 = pow(a,2.0);
   while (delta > toler) {
@@ -71,22 +109,26 @@ void run(double toler, double a)
       for (int j=1; j < M-1; j++) {
         phip[i][j] = (phi[i+1][j] + phi[i-1][j] +
                       phi[i][j+1] + phi[i][j-1])/4.0 +
-                      a2/(4.0 * epsilon0)*rhoa[i][j];
+                      a2*rhoa[i][j]/4.0;
       }
     }
     delta = mmax(phi, phip, M);
     swap(phi, phip, M);
 
   }
-  printf("iters %d", iter);
+  if(get_average_error(phi, analytical) > analytical_tol) {
+    printf("Failed to reach analytical solution, error %f\n", get_average_error(phi, analytical) );
+  }
+  return iter;
 
 }
 
 int main(int argc, char *argv[])
 {
+  int iter;
   clock_t start = clock();
-  run(1e-6, 0.01);
+  iter = run(1e-8, 1.0/(M-1));
   clock_t end = clock();
   double total = ((double)(end - start)) / CLOCKS_PER_SEC;
-  printf("Execution time: %f\n",total);
+  printf("Execution time %f, iters: %d \n",total, iter);
 }
diff --git a/poisson2d/optimized.f90 b/poisson2d/optimized.f90
deleted file mode 100644
index 3bc1045..0000000
--- a/poisson2d/optimized.f90
+++ /dev/null
@@ -1,55 +0,0 @@
-module rhofunc
-implicit none
-public
-integer, parameter :: dp=kind(0.d0)
-contains
-    pure real(dp) function rho(x,y)
-        real(dp), intent(in) :: x,y
-        if (x > 0.6_dp .and. x < 0.8_dp .and. y > 0.6_dp .and. y<0.8) then
-            rho = 1.0_dp
-        else if (x> 0.2_dp .and. x<0.4_dp .and. y>0.2_dp .and. y<0.4_dp) then
-            rho = -1.0_dp
-        else
-            rho = 0.0_dp
-        end if
-    end function
-
-end module
-
-program poisson
-use rhofunc, only: rho
-implicit none
-integer, parameter :: dp=kind(0.d0), M=300
-integer            :: i,j, iter
-real(dp),parameter :: epsilon0=8.85E-12_dp, target=1E-6_dp, a=0.01
-real(dp)           :: delta, b, e, phiprime(M,M), phi(M,M), a2, rhoarr(M,M)
-
-
-delta = 1.0_dp
-iter = 0
-call cpu_time(b)
-phiprime(:,:) = 0.0_dp
-phi(:,:) = 0.0_dp
-do i=1, M
-    do j=1, M
-        rhoarr(i,j) = rho(i*a,j*a)
-    end do
-end do
-
-do while (delta > target )
-    iter = iter + 1
-    a2 = a**2.0_dp
-    do i=2, M-1
-        do j=2, M-1
-            phiprime(i,j) = (phi(i+1,j) + phi(i-1,j) + phi(i,j+1) + phi(i,j-1))/4.0_dp &
-            + a2/4.0_dp/epsilon0*rhoarr(i,j)
-        end do
-    end do
-    delta = maxval(abs(phiprime - phi))
-    phi = phiprime 
-
-
-end do
-call cpu_time(e)
-print *, e-b, iter
-end program
\ No newline at end of file
diff --git a/poisson2d/poisson.py b/poisson2d/poisson.py
index e1659c1..1feaf19 100644
--- a/poisson2d/poisson.py
+++ b/poisson2d/poisson.py
@@ -5,49 +5,58 @@
 import matplotlib.pyplot as plt
 import datetime
 
-# Constants
-M = 100         # Grid squares on a side
-V = 1.0         # Voltage at top wall
-target = 1e-6   # Target accuracy
-a = 0.01
-epsilon0 = 8.85e-12
+
+M = 100      
+a = 1/(M-1)
+target = 1e-8
+analytical_tol = 5.0e-3
 
 # Create arrays to hold potential values
-phi = np.zeros([M+1,M+1],float)
-#phi[0,:] = V
-phiprime = np.zeros([M+1,M+1],float)
+phi = np.zeros([M,M],float)
+phiprime = np.zeros([M,M],float)
 
+#Analytical solution to the problem
+sol = np.zeros([M,M],float)
+begin = datetime.datetime.now()
+def realsol(x,y):
+    return y*(1-y)*x**3
+for i in range(M):
+    for j in range(M):
+        sol[i,j] = realsol(a*i,a*j)
+
+#Boundary conditions
+for i in range(M):
+    phi[M-1,i] = i*a*(1-i*a)
+    phiprime[M-1,i] = i*a*(1-i*a)
+
+#Source term
 def ro(x,y):
-    if x>0.6 and x<0.8 and y>0.6 and y<0.8:
-        return 1
-    elif x>0.2 and x<0.4 and y>0.2 and y<0.4:
-        return -1
-    else:
-        return 0
+    return -(6*x*y*(1-y) - 2*x**3)
 
 
 # Main loop
-begin = datetime.datetime.now()
 delta = 1.0
 while delta>target:
 
     # Calculate new values of the potential
     a2 = a**2
-    for i in range(1,M):
-        for j in range(1,M):
+    for i in range(1,M-1):
+        for j in range(1,M-1):
 
             phiprime[i,j] = (phi[i+1,j] + phi[i-1,j] \
                              + phi[i,j+1] + phi[i,j-1])/4 \
-                    + a2/4/epsilon0*ro(i*a,j*a)																												
+                    + a2*ro(i*a,j*a)/4																											
 
     # Calculate maximum difference from old values
     delta = np.max(abs(phi-phiprime))
 
     # Swap the two arrays around
     phi,phiprime = phiprime,phi
+maxsol = np.max(abs(sol))
+maxphi = np.max(abs(phi))
+sol = sol/maxsol
+phi = phi/maxphi
+assert np.sum(np.abs(sol-phi))/(M*M) < analytical_tol
 end = datetime.datetime.now()
 dif = end-begin
-print(dif.total_seconds())
-# Make a plot
-plt.imshow(phi,origin='lower')
-plt.savefig("purepython.png")
+print(dif.total_seconds())
\ No newline at end of file
diff --git a/poisson2d/poisson.pyx b/poisson2d/poisson.pyx
index 64d827f..d813d7c 100644
--- a/poisson2d/poisson.pyx
+++ b/poisson2d/poisson.pyx
@@ -9,59 +9,69 @@ import numpy as np
 cimport numpy as np
 import matplotlib.pyplot as plt
 cimport cython
+from libc.math cimport pow
 import time
 import datetime
 
 
-cdef int ro(double x,double y):
-    cdef double startone, endone, starttwo, endtwo
-    startone = 0.6
-    endone = 0.8
-    starttwo = 0.2
-    endtwo = 0.4
-    if x>startone and x<endone and y>startone and y<endone:
-        return 1
-    elif x>starttwo and x<endtwo and y>starttwo and y<endtwo:
-        return -1
-    else:
-        return 0
-
-
+cdef double ro(double x,double y):
+    cdef double res
+    res = -(6.0*x*y*(1.0-y)-2.0*pow(x,3.0))
+    return res
+cdef double boundary(double y):
+    return y*(1-y)
+cdef double analytical_solution(double x, double y):
+    return y*(1-y)*x**3
 @cython.cdivision(True)
 @cython.boundscheck(False)
 cdef void run_calculation():
     cdef int M, counter
-    cdef double V,target,a,epsilon0
-    M = 300         # Grid squares on a side
-    V = 1.0         # Voltage at top wall
-    target = 1e-6   # Target accuracy
-    a = 0.01
-    epsilon0 = 8.85e-12
+    cdef double V,target,a, analytical_tol
+    M = 100         # Grid squares on a side
+    target = 1e-8   # Target accuracy
+    a = 1.0/(M-1)
+    analytical_tol = 5.0e-3
 
     counter = 0
     cdef np.ndarray[np.double_t,ndim=2] phi = \
-        np.zeros([M+1,M+1],dtype=np.double)
+        np.zeros([M,M],dtype=np.double)
 
     cdef np.ndarray[np.double_t,ndim=2] phiprime = \
-        np.zeros([M+1,M+1],np.double)
+        np.zeros([M,M],np.double)
+    cdef np.ndarray[np.double_t,ndim=2] roa = \
+        np.zeros([M,M],np.double)
+    cdef np.ndarray[np.double_t,ndim=2] analytical_sol = \
+        np.zeros([M,M],np.double)
     cdef int i,j
     cdef double delta
     cdef double a2
     a2 = a**2
     delta = 1.0
+
+    for i in range(M):
+        for j in range(M):
+            roa[i,j] = ro(i*a, j*a)
+    for i in range(M):
+        phiprime[M-1,i] = boundary(i*a)
+        phi[M-1,i] = boundary(i*a)
+    for i in range(M):
+        for j in range(M):
+            analytical_sol[i,j] = analytical_solution(i*a,j*a)
+
     while delta>target:
         counter += 1
-        for i in range(1,M):
-            for j in range(1,M):
+        for i in range(1,M-1):
+            for j in range(1,M-1):
                 phiprime[i,j] = (phi[i+1,j] + phi[i-1,j] \
                                  + phi[i,j+1] + phi[i,j-1])/4 \
-                        + a2/(4*epsilon0)*ro(i*a,j*a)
+                        + a2*roa[i,j]/4.0
 
     # Calculate maximum difference from old values
-        delta = np.max(np.abs(phi-phiprime))
-
+        delta = np.max(abs(phi-phiprime))
     # Swap the two arrays around
         phi,phiprime = phiprime,phi
+    assert np.sum(np.abs(analytical_sol-phi))/(M*M) < analytical_tol
+    
     print(counter)
     
 begin = datetime.datetime.now()
diff --git a/poisson2d/poissonvectorized.py b/poisson2d/poissonvectorized.py
new file mode 100644
index 0000000..90233d5
--- /dev/null
+++ b/poisson2d/poissonvectorized.py
@@ -0,0 +1,64 @@
+#!/usr/bin/env python3
+# -*- coding: utf-8 -*-
+
+import numpy as np
+import matplotlib.pyplot as plt
+import datetime
+
+
+M = 300      
+a = 1/(M-1)
+target = 1e-8
+analytical_tol = 5.0e-3
+
+# Create arrays to hold potential values
+phi = np.zeros([M,M],float)
+phiprime = np.zeros([M,M],float)
+roa = np.zeros([M,M],float)
+#Analytical solution to the problem
+sol = np.zeros([M,M],float)
+begin = datetime.datetime.now()
+def realsol(x,y):
+    return y*(1-y)*x**3
+for i in range(M):
+    for j in range(M):
+        sol[i,j] = realsol(a*i,a*j)
+
+#Boundary conditions
+for i in range(M):
+    phi[M-1,i] = i*a*(1-i*a)
+    phiprime[M-1,i] = i*a*(1-i*a)
+
+#Source term
+def ro(x,y):
+    return -(6*x*y*(1-y) - 2*x**3)
+
+for i in range(M):
+    for j in range(M):
+        roa[i,j] = ro(i*a,j*a)
+
+# Main loop
+delta = 1.0
+iterations = 0
+while delta>target:
+
+    # Calculate new values of the potential
+    a2 = a**2
+    iterations += 1
+    phiprime[1:M-1,1:M-1] = (phi[2:M,1:M-1] + phi[0:M-2,1:M-1] \
+                     + phi[1:M-1,0:M-2] + phi[1:M-1,2:M])/4 \
+                    + a2*roa[1:M-1,1:M-1]/4																											
+
+    # Calculate maximum difference from old values
+    delta = np.max(abs(phi-phiprime))
+
+    # Swap the two arrays around
+    phi,phiprime = phiprime,phi
+maxsol = np.max(abs(sol))
+maxphi = np.max(abs(phi))
+sol = sol/maxsol
+phi = phi/maxphi
+assert np.sum(np.abs(sol-phi))/(M*M) < analytical_tol
+end = datetime.datetime.now()
+dif = end-begin
+print(dif.total_seconds(),iterations)
\ No newline at end of file
diff --git a/poisson2d/problemdesc.pdf b/poisson2d/problemdesc.pdf
new file mode 100755
index 0000000..4b67a46
Binary files /dev/null and b/poisson2d/problemdesc.pdf differ
diff --git a/poisson2d/problemdesc.tex b/poisson2d/problemdesc.tex
new file mode 100755
index 0000000..7120749
--- /dev/null
+++ b/poisson2d/problemdesc.tex
@@ -0,0 +1,38 @@
+\title{Problem: 2D Poisson equation}
+
+\date{\today}
+
+\documentclass[12pt]{article}
+\usepackage{amsmath}
+\usepackage[english]{babel}
+\begin{document}
+\maketitle
+\section{Theory}
+Partial differential equations are a mainstay of physics, and as such, many methods have been developed for solving them numerically. The problem to be solved in this idiom is obtaining a numerical solution to Poisson's equation
+\begin{align}
+\nabla ^2 u(x,y) = v(x,y)
+\end{align}
+by the finite-difference Jacobi relaxation method, i.e. by making a grid with some grid cell area $h^2$ and iteratively solving
+
+\begin{align}
+u_{i+1}(x,y) = &\frac{1}{4h^2}\bigg( u_i(x+h,y) + u_i(x-h,y)\\
+ &+ u_i(x,y-h) + u_i(x,y+h) - 4u_i(x,y)\bigg) - \frac{h^2}{4}v(x,y)
+\end{align}
+until some desired convergence. In other words, each grid point is updated by considering it and its neighbours in the previous iteration as well as the source term.
+\section{Specifics}
+\begin{itemize}
+\item The grid must be a unit square, i.e. the sides are of length 1. The grid spacing $h$ is $1/(M-1)$, if the number of grid cells is $M^2$.
+\item You must use a handwritten version of the algorithm above. Other than that, anything goes.
+\item The source term is $v(x,y) = 6xy(1-y)-2x^3$
+\item The initial guess for $u(x,y)$ should be zero everywhere except the boundary condition at $x=1$.
+\item The boundary conditions are $u(0,y) = 0, u(1,y) = y(1-y), u(x,0) = 0, u(x,1) = 0$.
+\item The code is convergent when the maximum grid point absolute difference between two iterations is $10^{-8}$. In other words, in pseudocode, when max(abs(current-previous))$<1e-8$.
+\item The exercise is considered complete when the convergent numerical solution has been checked against the analytical solution $u(x,y) = y(1-y)x^3$ and found to be within acceptable bounds. This is done by considering the average error per grid point, which must be below $5\cdot 10^{-3}$. Remember that the solution is unique only up to an additive constant, so you should normalize so that the highest absolute value of your solution (as well as the analytical solution) is 1!
+
+\end{itemize}
+
+Following these instructions, run the code for M=100, 200, 300. Your code should converge in 16054, 53666 and 106405 iterations, respectively; if it doesn't, check that your algorithm and initial conditions are correct. 
+\section{Notes}
+Be sure to be careful with the boundary conditions and take some care to set the array for $v(x,y)$ correctly. It is easy to make off-by-one errors that make convergence to the analytical solution impossible.
+\end{document}
+This is never printed
diff --git a/poisson2d/setup.py b/poisson2d/setup.py
deleted file mode 100644
index d732f30..0000000
--- a/poisson2d/setup.py
+++ /dev/null
@@ -1,16 +0,0 @@
-#!/usr/bin/env python3
-# -*- coding: utf-8 -*-
-"""
-Created on Thu Aug 31 16:09:12 2017
-
-@author: arc
-"""
-
-from distutils.core import setup
-from Cython.Build import cythonize
-import numpy
-setup(
-  name = 'yoo',
-  ext_modules = cythonize("poisson.pyx", annotate=True, language_level=3),
-  include_dirs=[numpy.get_include()],
-)
diff --git a/poisson2d/traditional1.f90 b/poisson2d/traditional1.f90
new file mode 100644
index 0000000..898d355
--- /dev/null
+++ b/poisson2d/traditional1.f90
@@ -0,0 +1,80 @@
+module naiveimpl
+implicit none
+public
+integer, parameter :: dp=kind(0.d0)
+contains
+
+pure real(dp) function rho(x,y)
+    real(dp), intent(in) :: x,y
+    rho = -(6.0_dp*x*y*(1.0_dp-y)-2.0_dp*x**3.0_dp)
+end function
+
+pure elemental real(dp) function boundary(y)
+    real(dp), intent(in) :: y
+
+    boundary = y*(1-y)
+end function
+
+pure real(dp) function analytical_sol(x,y)
+        real(dp), intent(in) :: x,y
+        analytical_sol = y*(1-y)*x**3
+end function
+
+integer function run_naive(M)
+    integer, intent(in):: M
+    integer, parameter :: dp=kind(0.d0)
+    real(dp)           :: phiprime(0:M-1,0:M-1), phi(0:M-1,0:M-1), a2, analytical(0:M-1,0:M-1), temp(0:M-1,0:M-1),delta,a, avgerror
+    real(dp),parameter :: target=1E-8_dp, analytical_tol = 5.0e-3
+    integer            :: i,j, iter
+    a = 1.0/(M-1)
+    delta = 1.0_dp
+    iter = 0
+    phi = 0
+    phiprime = 0
+    phiprime(M-1,:) = boundary([(i*a,i=0,M-1)])
+    phi(M-1,:) = boundary([(i*a,i=0,M-1)])
+    do i=0,M-1
+        do j=0, M-1
+            analytical(i,j) = analytical_sol(i*a,j*a)
+        end do
+    end do
+    do while (delta > target )
+        iter = iter + 1
+        a2 = a**2.0_dp
+        do i=1, M-2
+            do j=1, M-2
+                phiprime(i,j) = (phi(i+1,j) + phi(i-1,j) + phi(i,j+1) + phi(i,j-1))/4.0_dp &
+                + a2*rho(i*a,j*a)/4.0_dp
+            end do
+        end do
+        delta = maxval(abs(phiprime - phi))
+        temp = phi
+        phi = phiprime 
+        phiprime = temp
+    
+    end do
+    phi = phi/maxval(abs(phi))
+    analytical = analytical/maxval(abs(analytical))
+    avgerror = sum(abs(phi-analytical))/(M*M)
+    print *, avgerror
+    if(avgerror > analytical_tol) then
+        print *, "Failed to agree with analytical solution", avgerror
+    end if
+    run_naive = iter
+end function
+end module naiveimpl
+
+
+program poisson
+use naiveimpl, only: run_naive
+implicit none
+integer, parameter :: dp=kind(0.d0), M=100 
+integer            :: iter
+real(dp)           :: b,e
+
+
+call cpu_time(b)
+iter = run_naive(M)
+call cpu_time(e)
+print *, e-b, iter
+end program
\ No newline at end of file
diff --git a/poisson2d/traditional2.f90 b/poisson2d/traditional2.f90
new file mode 100644
index 0000000..3f63a9f
--- /dev/null
+++ b/poisson2d/traditional2.f90
@@ -0,0 +1,82 @@
+module traditionalopt
+implicit none
+public
+integer, parameter :: dp=kind(0.d0)
+contains
+pure real(dp) function rho(x,y)
+    real(dp), intent(in) :: x,y
+    rho = -(6.0_dp*x*y*(1.0_dp-y)-2.0_dp*x**3.0_dp)
+end function
+
+pure elemental real(dp) function boundary(y)
+    real(dp), intent(in) :: y
+
+    boundary = y*(1-y)
+end function
+
+pure real(dp) function analytical_sol(x,y)
+        real(dp), intent(in) :: x,y
+        analytical_sol = y*(1-y)*x**3
+end function
+
+integer function run_trad_optim(M)
+    integer            :: i,j, iter,M
+    real(dp),parameter ::  target=1E-8_dp ,analytical_tol = 5.0e-3
+    real(dp)           :: phiprime(0:M-1,0:M-1), phi(0:M-1,0:M-1), a2, analytical(0:M-1,0:M-1), source(0:M-1,0:M-1),delta,a,avgerror
+    a = 1.0/(M-1)
+    delta = 1.0_dp
+    iter = 0
+    phiprime(:,:) = 0.0_dp
+    phi(:,:) = 0.0_dp
+    phiprime(M-1,:) = boundary([(i*a,i=0,M-1)])
+    phi(M-1,:) = boundary([(i*a,i=0,M-1)])
+    do i=0,M-1
+        do j=0, M-1
+            analytical(i,j) = analytical_sol(i*a,j*a)
+        end do
+    end do
+
+    do i=0, M-1
+        do j=0, M-1
+            source(i,j) = rho(i*a,j*a)
+        end do
+    end do
+    
+
+    do while (delta > target )
+        iter = iter + 1
+        a2 = a**2.0_dp
+        do i=1, M-2
+            do j=1, M-2
+                phiprime(i,j) = (phi(i+1,j) + phi(i-1,j) + phi(i,j+1) + phi(i,j-1))/4.0_dp &
+                + a2*source(i,j)/4.0
+            end do
+        end do
+        delta = maxval(abs(phiprime - phi))
+        phi = phiprime 
+    end do
+    phi = phi/maxval(abs(phi))
+    analytical = analytical/maxval(abs(analytical))
+    avgerror = sum(abs(phi-analytical))/(M*M)
+    if(avgerror > analytical_tol) then
+        print *, "Failed to match analytical solution", avgerror
+    end if
+    run_trad_optim = iter
+    end function
+end module
+
+program poisson
+use traditionalopt, only: run_trad_optim
+implicit none
+integer, parameter :: dp=kind(0.d0), M=100
+integer            :: iter
+real(dp)           :: e,b
+
+
+
+
+call cpu_time(b)
+iter = run_trad_optim(M)
+call cpu_time(e)
+print *, e-b, iter
+end program
\ No newline at end of file