diff --git a/src/fhn_solver.jl b/src/fhn_solver.jl
index f22c072..2f8b7ec 100644
--- a/src/fhn_solver.jl
+++ b/src/fhn_solver.jl
@@ -24,48 +24,86 @@ using .Laplacian
     # Returns
     - `du`: calculated derivatives put back into the du array
 """
-function fhn!(du, u, p::CombinedPDEParams, t)
-    N = p.N
-    u_mat = reshape(u[1:N^2], N, N)
-    v_mat = reshape(u[N^2+1:end], N, N)
+# function fhn!(du, u, p::CombinedPDEParams, t)
+#     N = p.N
+#     u_mat = reshape(u[1:N^2], N, N)
+#     v_mat = reshape(u[N^2+1:end], N, N)
 
-    Δu = laplacian(u_mat, p.dx)
-    Δv = laplacian(v_mat, p.dx)
+#     Δu = laplacian(u_mat, p.dx)
+#     Δv = laplacian(v_mat, p.dx)
 
-    u_in = u_mat[2:end-1, 2:end-1]
-    v_in = v_mat[2:end-1, 2:end-1]
+#     u_in = u_mat[2:end-1, 2:end-1]
+#     v_in = v_mat[2:end-1, 2:end-1]
 
-    fu = p.Du * Δu .+ u_in .- u_in .^ 3 ./ 3 .- v_in
-    fv = p.Dv * Δv .+ p.ϵ * (u_in .+ p.a .- p.b .* v_in)
+#     fu = p.Du * Δu .+ u_in .- u_in .^ 3 ./ 3 .- v_in
+#     fv = p.Dv * Δv .+ p.ϵ * (u_in .+ p.a .- p.b .* v_in)
 
-    # Construct output with zero boundary (padding)
-    fu_full = zeros(N, N)
-    fv_full = zeros(N, N)
-    fu_full[2:end-1, 2:end-1] .= fu
-    fv_full[2:end-1, 2:end-1] .= fv
+#     # Construct output with zero boundary (padding)
+#     fu_full = zeros(N, N)
+#     fv_full = zeros(N, N)
+#     fu_full[2:end-1, 2:end-1] .= fu
+#     fv_full[2:end-1, 2:end-1] .= fv
 
-    du .= vcat(vec(fu_full), vec(fv_full))
-end
+#     du .= vcat(vec(fu_full), vec(fv_full))
+# end
 
-function step_fhn!(U, V, params_obs::Observable; dx=1)
+# function step_fhn!(U, V, params_obs::Observable; dx=1)
+#     p = params_obs[]
+
+#     # Flatten initial condition (activation u, recovery v)
+#     u0 = vec(U)
+#     v0 = vec(V)
+#     u_init = vcat(u0, v0)
+
+#     # Define one integration step using your fhn! function
+#     prob = ODEProblem((du, u, p, t) -> fhn!(du, u, p, t), u_init, (0.0, 0.1), p)
+#     sol = solve(prob, BS3(); save_start=false, saveat=10.0)
+
+#     # Extract solution and reshape
+#     u_new = reshape(sol[end][1:p.N^2], p.N, p.N)
+#     v_new = reshape(sol[end][p.N^2+1:end], p.N, p.N)
+
+#     # Update matrices in-place
+#     U .= u_new
+#     V .= v_new
+
+#     # Apply periodic boundary conditions
+#     U[1, :] .= U[end-1, :]
+#     U[end, :] .= U[2, :]
+#     U[:, 1] .= U[:, end-1]
+#     U[:, end] .= U[:, 2]
+
+#     V[1, :] .= V[end-1, :]
+#     V[end, :] .= V[2, :]
+#     V[:, 1] .= V[:, end-1]
+#     V[:, end] .= V[:, 2]
+
+#     return U
+# end
+function step_fhn!(U, V, params_obs::Observable; dx=1, dt=0.01)
     p = params_obs[]
 
-    # Flatten initial condition (activation u, recovery v)
-    u0 = vec(U)
-    v0 = vec(V)
-    u_init = vcat(u0, v0)
+    N = p.N
 
-    # Define one integration step using your fhn! function
-    prob = ODEProblem((du, u, p, t) -> fhn!(du, u, p, t), u_init, (0.0, 0.1), p)
-    sol = solve(prob, BS3(); save_start=false, saveat=10.0)
+    # Compute Laplacians on the interior
+    Δu = laplacian(U, dx)
+    Δv = laplacian(V, dx)
 
-    # Extract solution and reshape
-    u_new = reshape(sol[end][1:p.N^2], p.N, p.N)
-    v_new = reshape(sol[end][p.N^2+1:end], p.N, p.N)
+    # Extract interior
+    u_in = U[2:end-1, 2:end-1]
+    v_in = V[2:end-1, 2:end-1]
 
-    # Update matrices in-place
-    U .= u_new
-    V .= v_new
+    # Compute update using FHN reaction-diffusion terms
+    fu = p.Du .* Δu .+ u_in .- u_in .^ 3 ./ 3 .- v_in
+    fv = p.Dv .* Δv .+ p.ϵ .* (u_in .+ p.a .- p.b .* v_in)
+
+    # Euler update
+    u_new = u_in .+ dt .* fu
+    v_new = v_in .+ dt .* fv
+
+    # Write back interior updates
+    U[2:end-1, 2:end-1] .= u_new
+    V[2:end-1, 2:end-1] .= v_new
 
     # Apply periodic boundary conditions
     U[1, :] .= U[end-1, :]