diff --git a/poisson2d/README.md b/poisson2d/jacobi/README.md
similarity index 100%
rename from poisson2d/README.md
rename to poisson2d/jacobi/README.md
diff --git a/poisson2d/jacobi/legacy.f90 b/poisson2d/jacobi/legacy.f90
new file mode 100644
index 0000000..d12c31d
--- /dev/null
+++ b/poisson2d/jacobi/legacy.f90
@@ -0,0 +1,95 @@
+program poisson
+    !! This program solves the Poisson equation on the unit square with
+    !! homogeneous Dirichlet boundary conditions. The Laplace operator
+    !! is discretized with the standard second-order accurate central
+    !! finite difference scheme. The resulting discrete problem is solved
+    !! using the Jacobi iterative solver. The numerical implementation
+    !! relies on relatively standard Fortran constructs and is very similar
+    !! to what would be asked from students in a numerical analysis class.
+    implicit none
+    !----- Sets the double precision kind -----
+    integer, parameter :: precision = 15, range = 307
+    integer, parameter :: dp = selected_real_kind(precision, range)
+ 
+    !----- Physical parameters -----
+    integer , parameter :: m = 300
+    !! Number of grid points in each direction. 
+    real(dp), parameter :: dx = 1.0_dp / (m-1)
+    !! Uniform grid size in each direction.
+    real(dp) :: rhs(m, m)
+    !! Right-hand side of the Poisson equation.
+    real(dp) :: phi(m, m)
+    !! Solution of the Poisson equation.
+
+    !----- Jacobi solver -----
+    real(dp), parameter :: tolerance = 1e-6_dp
+    !! Tolerance for the iterative Jacobi solver.
+    real(dp) :: phi_prime(m, m)
+    !! Working array.
+    real(dp) :: residual
+    !! Residual is defined as the infinity-norm of the update, i.e. max | phi - phi' |
+    integer  :: iteration
+
+    !----- Miscellaneous -----
+    real(dp) :: start_time, end_time
+    integer             :: i, j
+    real(dp), parameter :: epsilon0 = 8.85E-12_dp
+
+    ! Initialize right-hand side.
+    do j = 1, m
+        do i = 1, m
+            rhs(i,j) = rho(i*dx,j*dx)
+        end do
+    end do
+
+    ! Tic.
+    call cpu_time(start_time)
+
+    !---------------------------------
+    !-----     JACOBI SOLVER     -----
+    !---------------------------------
+
+    ! Iteration counter.
+    iteration = 0
+    ! Working arrays.
+    phi = 0.0_dp ; phi_prime = 0.0_dp
+    ! Residual (ensuring at least one iteration is performed).
+    residual = 1.0_dp
+    ! Rescale rhs.
+    rhs = dx**2/(4*epsilon0) * rhs
+
+    ! Iterative solver.
+    do while (residual > tolerance)
+        ! Update iteration counter.
+        iteration = iteration + 1
+        ! Reset residual.
+        residual = 0.0_dp
+        ! Jacobi iteration.
+        do j = 2, m-1
+            do i = 2, m-1
+                ! Jacobi update.
+                phi_prime(i,j) = (phi(i+1,j)+phi(i-1,j)+phi(i,j+1)+phi(i,j-1))/4 + rhs(i, j)
+                ! On-the-fly residual computation.
+                residual = max(residual, abs(phi_prime(i,j) - phi(i,j)))
+            end do
+        end do
+        ! Updated solution.
+        phi = phi_prime 
+    end do
+   
+    ! Toc.
+    call cpu_time(end_time)
+
+    write(*, *) "Number of iterations  :", iteration
+    write(*, *) "Wall clock time (sec) :", end_time-start_time
+
+contains
+     pure real(dp) function rho(x, y)
+        real(dp), intent(in) :: x, y
+        if (all([x, y] > 0.6_dp .and. [x, y] < 0.8_dp)) then
+            rho = 1.0_dp
+        else
+            rho = merge(-1.0_dp, 0.0_dp, all([x, y]>0.2_dp .and. [x, y] < 0.4_dp))
+        endif
+   end function
+end program
diff --git a/poisson2d/naive.c b/poisson2d/jacobi/naive.c
similarity index 100%
rename from poisson2d/naive.c
rename to poisson2d/jacobi/naive.c
diff --git a/poisson2d/optimized.c b/poisson2d/jacobi/optimized.c
similarity index 100%
rename from poisson2d/optimized.c
rename to poisson2d/jacobi/optimized.c
diff --git a/poisson2d/poisson.py b/poisson2d/jacobi/poisson.py
similarity index 58%
rename from poisson2d/poisson.py
rename to poisson2d/jacobi/poisson.py
index e1659c1..6d5f051 100644
--- a/poisson2d/poisson.py
+++ b/poisson2d/jacobi/poisson.py
@@ -13,14 +13,15 @@
 epsilon0 = 8.85e-12
 
 # 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+1, M+1], float)
+# phi[0,:] = V
+phiprime = np.zeros([M+1, M+1], float)
 
-def ro(x,y):
-    if x>0.6 and x<0.8 and y>0.6 and y<0.8:
+
+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:
+    elif x > 0.2 and x < 0.4 and y > 0.2 and y < 0.4:
         return -1
     else:
         return 0
@@ -29,25 +30,25 @@ def ro(x,y):
 # Main loop
 begin = datetime.datetime.now()
 delta = 1.0
-while delta>target:
+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):
+        for j in range(1, M):
 
-            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)																												
+            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)
 
     # Calculate maximum difference from old values
     delta = np.max(abs(phi-phiprime))
 
     # Swap the two arrays around
-    phi,phiprime = phiprime,phi
+    phi, phiprime = phiprime, phi
 end = datetime.datetime.now()
 dif = end-begin
 print(dif.total_seconds())
 # Make a plot
-plt.imshow(phi,origin='lower')
+plt.imshow(phi, origin='lower')
 plt.savefig("purepython.png")
diff --git a/poisson2d/poisson.pyx b/poisson2d/jacobi/poisson.pyx
similarity index 100%
rename from poisson2d/poisson.pyx
rename to poisson2d/jacobi/poisson.pyx
diff --git a/poisson2d/setup.py b/poisson2d/jacobi/setup.py
similarity index 63%
rename from poisson2d/setup.py
rename to poisson2d/jacobi/setup.py
index d732f30..6b0104f 100644
--- a/poisson2d/setup.py
+++ b/poisson2d/jacobi/setup.py
@@ -10,7 +10,7 @@
 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()],
+    name='yoo',
+    ext_modules=cythonize("poisson.pyx", annotate=True, language_level=3),
+    include_dirs=[numpy.get_include()],
 )
diff --git a/poisson2d/naive.f90 b/poisson2d/naive.f90
deleted file mode 100644
index ecbc75a..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_dp) 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_dp
-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
diff --git a/poisson2d/optimized.f90 b/poisson2d/optimized.f90
deleted file mode 100644
index 38dddce..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_dp) 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_dp
-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