diff --git a/scripts/main.jl b/scripts/main.jl
index e5566c6..e38e373 100644
--- a/scripts/main.jl
+++ b/scripts/main.jl
@@ -1,6 +1,6 @@
 include("../src/gray_scott_solver.jl")
 include("../src/visualization.jl")
-include("../src/constants.jl")
+include("../src/utils/constants.jl")
 
 
 using Observables
diff --git a/scripts/run_simulation.jl b/scripts/run_simulation.jl
deleted file mode 100644
index 226ccf6..0000000
--- a/scripts/run_simulation.jl
+++ /dev/null
@@ -1,208 +0,0 @@
-using GLMakie, Observables
-include("../src/constants.jl")
-using .Constants
-include("../src/visualization.jl")
-using .Visualization
-
-# # Gray Scott Model Parameters
-# Parameters and initial conditions
-N = 128
-dx = 1.0
-# diffusion rates of substance 'u' and 'v'
-Du, Dv = Observable(0.16), Observable(0.08)
-# feed rate of 'u' and kill rate of 'v'
-F, k = Observable(0.060), Observable(0.062)
-dt = 1.0
-params_obs = Observable(GSParams(N, dx, Du[], Dv[], F[], k[]))
-
-# # GUI Parameters
-run_label = Observable{String}("Run")
-stepsize = Observable{Int}(30)
-
-
-
-function update_params!(params_obs::Observable, u, v, feed, kill)
-    old = params_obs[]
-    params_obs[] = GSParams(old.N, old.dx, u, v, feed, kill)
-end
-
-# Whenever a value gets changed via the textbox, update params object to reflect changes in simulation
-lift(Du, Dv, F, k) do u, v, F, k
-    update_params!(params_obs, u, v, F, k)
-end
-U = ones(N, N)
-V = zeros(N, N)
-center = N ÷ 2
-radius = 10
-# set a cube in the center with starting concentrations for 'u' and 'v'
-U[center-radius:center+radius, center-radius:center+radius] .= 0.50
-V[center-radius:center+radius, center-radius:center+radius] .= 0.25
-
-# Observable holding current U for heatmap
-heat_obs = Observable(U)
-
-
-function laplacian5(f)
-    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
-end
-
-function step!(U, V, params_obs::Observable)
-    lap_u = laplacian5(U)
-    lap_v = laplacian5(V)
-    diff_u = params_obs[].Du
-    diff_v = params_obs[].Dv
-    feed_u = params_obs[].F
-    kill_v = params_obs[].k
-    u = U[2:end-1, 2:end-1]
-    v = V[2:end-1, 2:end-1]
-
-    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)
-
-    # Update with new values
-    U[2:end-1, 2:end-1] .= u_new
-    V[2:end-1, 2:end-1] .= v_new
-
-    # 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]
-
-    # Update heatmap observable
-    return U
-end
-function multi_step!(state, n_steps)
-    for _ in 1:n_steps
-        heat_obs[] = step!(state[1], state[2], params_obs)
-    end
-end
-
-# Build GUI
-fig = Figure(size=(600, 600))
-
-
-gh = GridLayout(fig[1, 1])
-ax = Axis(gh[1, 1])
-
-hm = Makie.heatmap!(ax, heat_obs, colormap=:viridis)
-Makie.deactivate_interaction!(ax, :rectanglezoom)
-ax.aspect = DataAspect()
-
-spoint = select_point(ax.scene)
-
-function coord_to_index(x, y, N)
-    ix = clamp(round(Int, x), 1, N)
-    iy = clamp(round(Int, y), 1, N)
-    return ix, iy
-end
-r = 5
-
-on(spoint) do pt
-    if pt === nothing
-        return
-    end
-    x, y = pt
-    i, j = coord_to_index(x, y, N)
-
-    # get corners of square that will get filled with concentration
-    imin = max(i - r, 1)
-    imax = min(i + r, N)
-    jmin = max(j - r, 1)
-    jmax = min(j + r, N)
-
-    # set disbalanced concentration of U and V
-    U[imin:imax, jmin:jmax] .= 0.5
-    V[imin:imax, jmin:jmax] .= 0.25
-
-    heat_obs[] = copy(U)
-end
-
-# # Controls
-fig[2, 1] = buttongrid = GridLayout(ax.scene, tellwidth=false)
-btn_step = Button(buttongrid[1, 1], width=50, label="Step")
-btn_start = Button(buttongrid[1, 2], width=50, label=run_label)
-btn_reset = Button(buttongrid[1, 3], width=50, label="Reset")
-slidergrid = SliderGrid(fig[3, 1], (label="Speed", range=1:1:100, format="{}x", width=350, startvalue=stepsize[]))
-
-speed_slider = slidergrid.sliders[1].value
-on(speed_slider) do s
-    try
-        stepsize[] = s
-        println("Changed stepsize to $s")
-    catch
-        @warn "Invalid input for $s"
-    end
-end
-
-gh[1, 2] = textboxgrid = GridLayout(ax.scene, tellwidth=false)
-rowsize!(gh, 1, Relative(1.0))
-function param_box!(row, labeltxt, observable)
-    Label(textboxgrid[row, 1], labeltxt)
-    box = Textbox(textboxgrid[row, 2], validator=Float64, width=50, placeholder=labeltxt, stored_string="$(observable[])")
-    on(box.stored_string) do s
-        try
-            observable[] = parse(Float64, s)
-            println("changed $labeltxt to $s")
-        catch
-            @warn "Invalid input for $labeltxt: $s"
-        end
-    end
-end
-param_box!(1, "Du", Du)
-param_box!(2, "Dv", Dv)
-param_box!(3, "Feed", F)
-param_box!(4, "kill", k)
-
-
-# Timer and state for animation
-running = Observable(false)
-function reset!(U, V, heat_obs)
-    U .= 1.0
-    V .= 0.0
-    center = size(U, 1) ÷ 2
-    radius = 10
-    U[center-radius:center+radius, center-radius:center+radius] .= 0.50
-    V[center-radius:center+radius, center-radius:center+radius] .= 0.25
-    heat_obs[] = copy(U)
-end
-
-on(running) do r
-    run_label[] = r ? "Pause" : "Run"
-end
-
-
-# Button Listeners
-on(btn_step.clicks) do _
-    multi_step!((U, V), stepsize[])
-end
-
-on(btn_start.clicks) do _
-    running[] = !running[]
-end
-
-on(btn_start.clicks) do _
-    @async while running[]
-        multi_step!((U, V), stepsize[])
-        sleep(0.05)  # ~20 FPS
-    end
-end
-
-on(btn_reset.clicks) do _
-    running[] = false
-    reset!(U, V, heat_obs)
-
-end
-display(fig)
diff --git a/src/gray_scott_solver.jl b/src/gray_scott_solver.jl
index ca1840b..3668d26 100644
--- a/src/gray_scott_solver.jl
+++ b/src/gray_scott_solver.jl
@@ -1,19 +1,10 @@
 module GrayScottSolver
 using Observables
-include("constants.jl")
+include("utils/constants.jl")
+include("utils/laplacian.jl")
 using .Constants
-export laplacian5, step!, multi_step!
-
-
-function laplacian5(f, dx)
-    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
-end
+using .Laplacian
+export step!, multi_step!
 
 function step!(U, V, params_obs::Observable; dx=1)
     lap_u = laplacian5(U, dx)
diff --git a/src/solver.jl b/src/solver.jl
deleted file mode 100644
index 3910c8f..0000000
--- a/src/solver.jl
+++ /dev/null
@@ -1,154 +0,0 @@
-using DifferentialEquations
-using Random
-include("constants.jl")
-using .Constants
-
-"""
-    fhn(du, u, p:FHNParams, t:)
-
-    Implements the spatial dynamics of FitzHugh-Nagumo (fhn). Designed to be 
-    within a larger numerical solver of partial differential equations.
-
-    # Arguments
-    - `du`: output argument which stores the calculated derivatives 
-    - `u`: input vector containing the current state of the system at time t
-    - `p`: holds all the fixed parameters of the FHN model
-    - `t`: current time 
-
-    # Returns
-    - `du`: calculated derivatives put back into the du array
-"""
-function fhn!(du, u, p::FHNParams, t=0)
-    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
-
-    Δ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)
-
-    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)
-
-    du .= vcat(vec(fu), vec(fv))
-end
-
-"""
-    run_simulation(tspan::Tuple{Float64,Float64}, N::Int)
-
-    solving the ODE and modelling it after FHN
-
-    # Arguments
-    - `tspan`: tuple of two Float64's representing start and end times for simulation
-    - `N`: size of the N×N grid
-
-    # Returns
-    - `sol`: solved differential equation (ODE)
-"""
-function run_simulation(tspan::Tuple{Float64,Float64}, N::Int)
-    # Turing-spot parameters
-    p = FHNParams(N=N, dx=1.0, Du=1e-5, Dv=1e-3, ϵ=0.01, a=0.1, b=0.5)
-
-    # Initial conditions (random noise)
-    Random.seed!(1234)
-    # Create two vectors with length N*N with numbers between 0.1 and 0.11
-    u0 = vec(0.1 .+ 0.01 .* rand(N, N))
-    v0 = vec(0.1 .+ 0.01 .* rand(N, N))
-
-    y0 = vcat(u0, v0)
-
-    prob = ODEProblem(fhn!, y0, tspan, p)
-    sol = solve(prob, BS3(), saveat=1.0)  # You can try `Rosenbrock23()` too
-
-    return sol
-end
-
-# Laplacian with 5-point stencil
-function laplacianA(A::Matrix{Float64}, N::Int, dx::Float64)
-    dx2 = dx^2
-    lap = zeros(N, N)
-    lap[2:end-1, 2:end-1] .= (
-        A[1:end-2, 2:end-1] .+ A[3:end, 2:end-1]  # up/down
-        .+
-        A[2:end-1, 1:end-2] .+ A[2:end-1, 3:end]  # left/right
-        .-
-        4 .* A[2:end-1, 2:end-1]
-    ) ./ dx2
-    return lap
-end
-
-# Gray-Scott model: reaction + diffusion
-function grayscott!(du, u, p::GSParams, t=0)
-    N = p.N
-    U = reshape(u[1:N^2], N, N)
-    V = reshape(u[N^2+1:end], N, N)
-
-    ΔU = laplacianA(U, N, p.dx)
-    ΔV = laplacianA(V, N, p.dx)
-
-    UVV = U .* V .* V
-
-    dU = p.Du * ΔU .- UVV .+ p.F .* (1 .- U)
-    dV = p.Dv * ΔV .+ UVV .- (p.F .+ p.k) .* V
-
-    du .= vcat(vec(dU), vec(dV))
-end
-
-# Run simulation
-function run_simulationG(tspan::Tuple{Float64,Float64}, N::Int, GSP::GSParams)
-    # Replace this in run_simulation
-    p = GSP
-
-    # Initial conditions: U = 1, V = 0 everywhere
-    U = ones(N, N)
-    V = zeros(N, N)
-
-    # Seed a square in the center
-    center = N ÷ 2
-    radius = 10
-    U[center-radius:center+radius, center-radius:center+radius] .= 0.50
-    V[center-radius:center+radius, center-radius:center+radius] .= 0.25
-
-    y0 = vcat(vec(U), vec(V))
-
-    prob = ODEProblem(grayscott!, y0, tspan, p)
-    sol = solve(prob, BS3(), saveat=10.0)  # or Rosenbrock23()
-
-    return sol
-end
-
-function run_simulationG_no_ode(tspan::Tuple{Float64,Float64}, N::Int, GSP::GSParams)
-    p = GSP
-    dx2 = p.dx^2
-    dt = 1.0  # fixed timestep
-
-    # Initialize U and V
-    U = ones(N, N)
-    V = zeros(N, N)
-
-    # Seed pattern in the center
-    center = N ÷ 2
-    radius = 10
-    U[center-radius:center+radius, center-radius:center+radius] .= 0.50
-    V[center-radius:center+radius, center-radius:center+radius] .= 0.25
-
-    # Store history for visualization (optional)
-    snapshots = []
-
-    for i in 1:tspan[2]
-        print("Running $i/$(tspan[2])\r")
-
-        lap_U = laplacianA(U, N, p.dx)
-        lap_V = laplacianA(V, N, p.dx)
-
-        UVV = U .* V .* V
-        dU = p.Du * lap_U .- UVV .+ p.F * (1 .- U)
-        dV = p.Dv * lap_V .+ UVV .- (p.F + p.k) .* V
-
-        U .+= dU * dt
-        V .+= dV * dt
-
-        if i % 100 == 0  # save every 100 steps
-            push!(snapshots, (copy(U), copy(V)))
-        end
-    end
-    return snapshots
-end
diff --git a/src/constants.jl b/src/utils/constants.jl
similarity index 100%
rename from src/constants.jl
rename to src/utils/constants.jl
diff --git a/src/laplacian.jl b/src/utils/laplacian.jl
similarity index 71%
rename from src/laplacian.jl
rename to src/utils/laplacian.jl
index f9f4965..e04d5a2 100644
--- a/src/laplacian.jl
+++ b/src/utils/laplacian.jl
@@ -1,3 +1,5 @@
+module Laplacian
+
 """
     laplacian(U::Matrix{Float64}, N::Int, h::Float64)
 
@@ -18,3 +20,17 @@ function laplacian(U::Matrix{Float64}, N::Int, h::Float64)
              circshift(U, (0, -1)) .+ circshift(U, (0, 1)) .- 4 .* U
     return vec(padded) ./ h^2
 end
+
+function laplacian5(f, dx)
+    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
+end
+
+export laplacian, laplacian5
+
+end # module Laplacian
\ No newline at end of file
diff --git a/src/visualization.jl b/src/visualization.jl
index 7d1a066..d2f8343 100644
--- a/src/visualization.jl
+++ b/src/visualization.jl
@@ -38,11 +38,13 @@ function step_through_solution(sol, N::Int)
 
     display(fig)
 end
+
 function coord_to_index(x, y, N)
     ix = clamp(round(Int, x), 1, N)
     iy = clamp(round(Int, y), 1, N)
     return ix, iy
 end
+
 function reset!(U, V, heat_obs)
     U .= 1.0
     V .= 0.0
@@ -52,6 +54,7 @@ function reset!(U, V, heat_obs)
     V[center-radius:center+radius, center-radius:center+radius] .= 0.25
     heat_obs[] = copy(U)
 end
+
 function param_box!(grid, row, labeltxt, observable::Observable)
     Label(grid[row, 1], labeltxt)
     box = Textbox(grid[row, 2], validator=Float64, width=50, placeholder=labeltxt, stored_string="$(observable[])")