diff --git a/scripts/main.jl b/scripts/main.jl
index dd0a440..b9f5213 100644
--- a/scripts/main.jl
+++ b/scripts/main.jl
@@ -1,6 +1,6 @@
 include("../src/utils/constants.jl")
-include("../src/fhn_solver.jl")
-#include("../src/gray_scott_solver.jl")
+#include("../src/fhn_solver.jl")
+include("../src/gray_scott_solver.jl")
 include("../src/visualization.jl")
 
 using Observables
@@ -49,12 +49,13 @@ function setup_param_binding!(model_type, gray_params, fhn_params, params_obs)
 end
 
 lift(Du, Dv, ϵ, a, b) do d_u, d_v, eps, aa, bb
-    params_obs[] = Constants.FHNParams(N=N, dx=dx, Du=d_u, Dv=d_v, ϵ=eps, a=aa, b=bb)
+    params_obs[] = FHNParams(N=N, dx=dx, Du=d_u, Dv=d_v, ϵ=eps, a=aa, b=bb)
 end
+"""
 
 U = ones(N, N)
 V = zeros(N, N)
 heat_obs = Observable(U)
 
-fig = Visualization.build_ui(U, V, param_observables, params_obs, heat_obs)
+fig = build_ui(U, V, param_observables, params_obs, heat_obs)
 display(fig)
diff --git a/src/fhn_solver.jl b/src/fhn_solver.jl
index 2dbc3cd..cd412e7 100644
--- a/src/fhn_solver.jl
+++ b/src/fhn_solver.jl
@@ -23,16 +23,26 @@ using .Laplacian
     - `du`: calculated derivatives put back into the du array
 """
 function fhn!(du, u, p::FHNParams, t)
-    u_mat = reshape(u[1:p.N^2], p.N, p.N) # activation variable
-    v_mat = reshape(u[p.N^2+1:end], p.N, p.N) # deactivation variable
+    N = p.N
+    u_mat = reshape(u[1:N^2], N, N)
+    v_mat = reshape(u[N^2+1:end], N, N)
 
-    Δu = reshape(laplacian(u_mat, p.N, p.dx), p.N, p.N)
-    Δv = reshape(laplacian(v_mat, p.N, p.dx), p.N, p.N)
+    Δu = laplacian(u_mat, p.dx)
+    Δv = laplacian(v_mat, p.dx)
 
-    fu = p.Du * Δu .+ u_mat .- u_mat .^ 3 ./ 3 .- v_mat
-    fv = p.Dv * Δv .+ p.ϵ * (u_mat .+ p.a .- p.b .* v_mat)
+    u_in = u_mat[2:end-1, 2:end-1]
+    v_in = v_mat[2:end-1, 2:end-1]
 
-    du .= vcat(vec(fu), vec(fv))
+    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
+
+    du .= vcat(vec(fu_full), vec(fv_full))
 end
 
 function step_fhn!(U, V, params_obs::Observable; dx=1)
diff --git a/src/gray_scott_solver.jl b/src/gray_scott_solver.jl
index 8fcb29d..ffc7161 100644
--- a/src/gray_scott_solver.jl
+++ b/src/gray_scott_solver.jl
@@ -7,24 +7,29 @@ using .Constants
 using .Laplacian
 
 function step_gray_scott!(U, V, params_obs::Observable; dx=1)
-    lap_u = laplacian5(U, dx)
-    lap_v = laplacian5(V, dx)
-    diff_u = params_obs[].Du
-    diff_v = params_obs[].Dv
-    feed_u = params_obs[].F
-    kill_v = params_obs[].k
+    # Extract parameters
+    p = params_obs[]
+    Du, Dv = p.Du, p.Dv
+    F, k = p.F, p.k
+
+    # Compute Laplacians on the interior
+    lap_u = laplacian(U, dx)  # shape (N-2, N-2)
+    lap_v = laplacian(V, dx)
+
+    # Extract interior values
     u = U[2:end-1, 2:end-1]
     v = V[2:end-1, 2:end-1]
 
+    # Gray-Scott update equations (Euler-style)
     uvv = u .* v .* v
-    u_new = u .+ (diff_u .* lap_u .- uvv .+ feed_u .* (1 .- u))
-    v_new = v .+ (diff_v .* lap_v .+ uvv .- (feed_u .+ kill_v) .* v)
+    u_new = u .+ (Du .* lap_u .- uvv .+ F .* (1 .- u))
+    v_new = v .+ (Dv .* lap_v .+ uvv .- (F .+ k) .* v)
 
-    # Update with new values
+    # Update interior
     U[2:end-1, 2:end-1] .= u_new
     V[2:end-1, 2:end-1] .= v_new
 
-    # Periodic boundary conditions
+    # Apply periodic boundary conditions
     U[1, :] .= U[end-1, :]
     U[end, :] .= U[2, :]
     U[:, 1] .= U[:, end-1]
@@ -35,8 +40,7 @@ function step_gray_scott!(U, V, params_obs::Observable; dx=1)
     V[:, 1] .= V[:, end-1]
     V[:, end] .= V[:, 2]
 
-    # Update heatmap observable
-    return U
+    return U  # for visualization
 end
 function multi_step!(state, n_steps, heat_obs::Observable, params_obs::Observable; dx=1)
     for _ in 1:n_steps
diff --git a/src/utils/laplacian.jl b/src/utils/laplacian.jl
index a9cebdb..814573a 100644
--- a/src/utils/laplacian.jl
+++ b/src/utils/laplacian.jl
@@ -1,35 +1,14 @@
 module Laplacian
 
-"""
-    laplacian(U::Matrix{Float64}, N::Int, h::Float64)
-
-    Computes the discrete 2D Laplacian of a matrix `U` using a 5-point stencil
-    and circular boundary conditions.
-
-    # Arguments
-    - `U::Matrix{Float64}`: The input 2D matrix representing the field or image.
-    - `N::Int`: Integer
-    - `h::Float64`: The spatial step size or grid spacing between points in the discretization.
-
-    # Returns
-    - `Vector{Float64}`: A flattened (vectorized) representation of the approximated Laplacian values for each element in `U`. The boundary conditions are handled circularly.
-"""
-function laplacian(U::Matrix{Float64}, N::Int, h::Float64)
-    Δ = zeros(N, N)
-    for i in 2:N-1, j in 2:N-1
-        Δ[i, j] = (U[i+1, j] + U[i-1, j] + U[i, j+1] + U[i, j-1] - 4 * U[i, j]) / h^2
-    end
-    return vec(Δ)
-end
-
-function laplacian5(f, dx)
+function laplacian(U::AbstractMatrix{<:Real}, dx::Real)
     h2 = dx^2
-    left = f[2:end-1, 1:end-2]
-    right = f[2:end-1, 3:end]
-    down = f[3:end, 2:end-1]
-    up = f[1:end-2, 2:end-1]
-    center = f[2:end-1, 2:end-1]
-    return (left .+ right .+ down .+ up .- 4 .* center) ./ h2
+    center = U[2:end-1, 2:end-1]
+    up     = U[1:end-2, 2:end-1]
+    down   = U[3:end,   2:end-1]
+    left   = U[2:end-1, 1:end-2]
+    right  = U[2:end-1, 3:end]
+
+    return (up .+ down .+ left .+ right .- 4 .* center) ./ h2
 end
 
 export laplacian, laplacian5
diff --git a/src/visualization.jl b/src/visualization.jl
index 930f288..751764a 100644
--- a/src/visualization.jl
+++ b/src/visualization.jl
@@ -1,6 +1,6 @@
 module Visualization
-#include("gray_scott_solver.jl")
-include("fhn_solver.jl")
+include("gray_scott_solver.jl")
+#include("fhn_solver.jl")
 
 using GLMakie, Observables, Makie