Manifest.toml

@@ -2,7 +2,7 @@
 
 julia_version = "1.11.4"
 manifest_format = "2.0"
-project_hash = "873df6eea0990fe68ca1cd8039db95353f9b799c"
+project_hash = "0d147eaf78630b05e95b5317275a880443440272"
 
 [[deps.ADTypes]]
 git-tree-sha1 = "e2478490447631aedba0823d4d7a80b2cc8cdb32"

Project.toml

@@ -1,3 +1,4 @@
 [deps]
 DifferentialEquations = "0c46a032-eb83-5123-abaf-570d42b7fbaa"
 GLMakie = "e9467ef8-e4e7-5192-8a1a-b1aee30e663a"
+Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"

src/initial_conditions.jl

@@ -1,3 +1,5 @@
+import Random
+
 function initial_conditions(N)
     Random.seed!(1234)
     u = 0.1 .+ 0.01 .* rand(N, N)

src/laplacian.jl

@@ -1,5 +1,4 @@
-function laplacian(u, N, h)
+function laplacian(U::Matrix{Float64}, N::Int, h::Float64)
-    U = reshape(u, N, N)
     padded = circshift(U, (-1, 0)) .+ circshift(U, (1, 0)) .+
              circshift(U, (0, -1)) .+ circshift(U, (0, 1)) .- 4 .* U
     return vec(padded) ./ h^2

src/solver.jl

@@ -1,34 +1,35 @@
+using DifferentialEquations
+
+
 function run_simulation(N::Int, tspan::Tuple{Float64,Float64})
-    # ✅ Turing pattern parameters — promotes stable spots
+    # Turing-spot parameters
     a = 0.1
     b = 0.5
     ϵ = 0.01
-    Du = 1e-5   # activator diffusion
+    Du = 1e-5
-    Dv = 1e-3   # inhibitor diffusion (>> Du)
+    Dv = 1e-3
 
-    # Laplacian setup (assumes ∆x = 1)
+    dx = 1.0  # grid spacing
-    dx = 1.0
-    L = laplacian_matrix(N, dx)
 
-    # ✅ Local random noise around a constant for initial conditions
+    # Initial conditions (random noise)
-    u0, v0 = initial_conditions(N)  # Should include a seeded RNG
+    u0, v0 = initial_conditions(N)
     y0 = vcat(u0, v0)
 
-    # ODE system for Reaction-Diffusion with FitzHugh-Nagumo dynamics
     function f!(du, u, p, t)
-        U = u[1:N^2]
+        u_mat = reshape(u[1:N^2], N, N)
-        V = u[N^2+1:end]
+        v_mat = reshape(u[N^2+1:end], N, N)
 
-        ΔU = Du .* (L * U)
+        Δu = reshape(laplacian(u_mat, N, dx), N, N)
-        ΔV = Dv .* (L * V)
+        Δv = reshape(laplacian(v_mat, N, dx), N, N)
 
-        du[1:N^2] = ΔU .+ (U .- (U .^ 3) ./ 3 .- V) ./ ϵ
+        fu = Du * Δu .+ u_mat .- u_mat .^ 3 ./ 3 .- v_mat
-        du[N^2+1:end] = ΔV .+ ϵ .* (U .+ a .- b .* V)
+        fv = Dv * Δv .+ ϵ * (u_mat .+ a .- b .* v_mat)
+
+        du .= vcat(vec(fu), vec(fv))
     end
 
-    # ✅ Use stiff solver for stability
     prob = ODEProblem(f!, y0, tspan)
-    sol = solve(prob, Rosenbrock23(), saveat=1.0, reltol=1e-6, abstol=1e-6)
+    sol = solve(prob, BS3(), saveat=1.0)  # You can try `Rosenbrock23()` too
 
     return sol
 end