module FHNSolver

include("utils/laplacian.jl")

using DifferentialEquations
using Random
using Observables

using .Laplacian

function step_fhn!(U, V, params_obs::Observable; dx=1, dt=0.01)
    p = params_obs[]

    N = p.N

    # Compute Laplacians on the interior
    Δu = laplacian(U, dx)
    Δv = laplacian(V, dx)

    # Extract interior
    u_in = U[2:end-1, 2:end-1]
    v_in = V[2:end-1, 2:end-1]

    # 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, :]
    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

end # Module end