diff --git a/.gitignore b/.gitignore
index c487889..7a4896f 100644
--- a/.gitignore
+++ b/.gitignore
@@ -24,3 +24,4 @@ docs/site/
 # environment.
 Manifest.toml
 
+gif
diff --git a/Manifest.toml b/Manifest.toml
index 23b9f57..eceefdd 100644
--- a/Manifest.toml
+++ b/Manifest.toml
@@ -1,8 +1,8 @@
 # This file is machine-generated - editing it directly is not advised
 
-julia_version = "1.11.5"
+julia_version = "1.11.4"
 manifest_format = "2.0"
-project_hash = "0d147eaf78630b05e95b5317275a880443440272"
+project_hash = "b5f5f0b50b1e0b7dd05bb93e0b1bb03fea235d53"
 
 [[deps.ADTypes]]
 git-tree-sha1 = "e2478490447631aedba0823d4d7a80b2cc8cdb32"
@@ -1639,7 +1639,7 @@ version = "3.2.4+0"
 [[deps.OpenLibm_jll]]
 deps = ["Artifacts", "Libdl"]
 uuid = "05823500-19ac-5b8b-9628-191a04bc5112"
-version = "0.8.5+0"
+version = "0.8.1+4"
 
 [[deps.OpenSSL_jll]]
 deps = ["Artifacts", "JLLWrappers", "Libdl"]
diff --git a/Project.toml b/Project.toml
index 19c0af6..0672cef 100644
--- a/Project.toml
+++ b/Project.toml
@@ -1,4 +1,6 @@
 [deps]
 DifferentialEquations = "0c46a032-eb83-5123-abaf-570d42b7fbaa"
 GLMakie = "e9467ef8-e4e7-5192-8a1a-b1aee30e663a"
+Makie = "ee78f7c6-11fb-53f2-987a-cfe4a2b5a57a"
+Observables = "510215fc-4207-5dde-b226-833fc4488ee2"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
diff --git a/scripts/main.jl b/scripts/main.jl
new file mode 100644
index 0000000..e38e373
--- /dev/null
+++ b/scripts/main.jl
@@ -0,0 +1,28 @@
+include("../src/gray_scott_solver.jl")
+include("../src/visualization.jl")
+include("../src/utils/constants.jl")
+
+
+using Observables
+using GLMakie
+using .GrayScottSolver
+using .Visualization
+using .Constants
+
+const N = 128
+const dx = 1.0
+Du, Dv = Observable(0.16), Observable(0.08)
+F, k = Observable(0.060), Observable(0.062)
+
+params_obs = Observable(GSParams(N, dx, Du[], Dv[], F[], k[]))
+lift(Du, Dv, F, k) do u, v, f, ki
+    params_obs[] = GSParams(N, dx, u, v, f, ki)
+end
+
+
+U = ones(N, N)
+V = zeros(N, N)
+heat_obs = Observable(U)
+
+fig = build_ui(U, V, Du, Dv, F, k, params_obs, heat_obs)
+display(fig)
diff --git a/scripts/run_simulation.jl b/scripts/run_simulation.jl
deleted file mode 100644
index 3e43bb4..0000000
--- a/scripts/run_simulation.jl
+++ /dev/null
@@ -1,13 +0,0 @@
-include("../src/AnimalFurFHN.jl")  # this loads the module code
-using .AnimalFurFHN
-# include optional visualizer only if needed:
-include("../src/visualization.jl")
-using .Visualization
-
-N = 100
-tspan = (0.0, 1000.0)
-
-sol = AnimalFurFHN.run_simulation(tspan, N)
-
-
-Visualization.step_through_solution(sol, N)
\ No newline at end of file
diff --git a/src/AnimalFurFHN.jl b/src/AnimalFurFHN.jl
deleted file mode 100644
index 060abd7..0000000
--- a/src/AnimalFurFHN.jl
+++ /dev/null
@@ -1,9 +0,0 @@
-# src/AnimalFurFHN.jl
-module AnimalFurFHN
-
-include("constants.jl")
-include("laplacian.jl")
-include("solver.jl")
-
-export run_simulation   # Make sure this is here!
-end
diff --git a/src/gray_scott_solver.jl b/src/gray_scott_solver.jl
new file mode 100644
index 0000000..3668d26
--- /dev/null
+++ b/src/gray_scott_solver.jl
@@ -0,0 +1,46 @@
+module GrayScottSolver
+using Observables
+include("utils/constants.jl")
+include("utils/laplacian.jl")
+using .Constants
+using .Laplacian
+export step!, multi_step!
+
+function step!(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
+    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, heat_obs::Observable, params_obs::Observable; dx=1)
+    for _ in 1:n_steps
+        heat_obs[] = step!(state[1], state[2], params_obs; dx=1)
+    end
+end
+end
\ No newline at end of file
diff --git a/src/solver.jl b/src/solver.jl
deleted file mode 100644
index 8999319..0000000
--- a/src/solver.jl
+++ /dev/null
@@ -1,61 +0,0 @@
-using DifferentialEquations
-using Random
-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
diff --git a/src/constants.jl b/src/utils/constants.jl
similarity index 65%
rename from src/constants.jl
rename to src/utils/constants.jl
index d1734e3..c06f46c 100644
--- a/src/constants.jl
+++ b/src/utils/constants.jl
@@ -15,6 +15,16 @@ struct FHNParams
         new(N, dx, Du, Dv, ϵ, a, b)
 end
 
-export FHNParams
+struct GSParams
+    N::Int       # grid size
+    dx::Float64  # grid spacing
+    Du::Float64  # diffusion rate U
+    Dv::Float64  # diffusion rate V
+    F::Float64   # feed rate
+    k::Float64   # kill rate
+
+end
+
+export FHNParams, GSParams
 
 end # module Constants
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 7f1514a..d2f8343 100644
--- a/src/visualization.jl
+++ b/src/visualization.jl
@@ -1,6 +1,7 @@
 module Visualization
-
-using GLMakie
+include("gray_scott_solver.jl")
+using GLMakie, Observables, Makie
+using .GrayScottSolver
 """
     step_through_solution(sol::SolutionType, N::Int)
 
@@ -18,18 +19,146 @@ function step_through_solution(sol, N::Int)
     ax = Axis(fig[1, 1])
     slider = Slider(fig[2, 1], range=1:length(sol), startvalue=1)
 
+    textbox = Textbox(fig[1, 2], placeholder="Feed Rate", validator=Float64)
+    F = Observable(0.060)
     # Initialize heatmap with first time step
     u0 = reshape(sol[1][1:N^2], N, N)
     heat_obs = Observable(u0)
     hmap = heatmap!(ax, heat_obs, colormap=:magma)
 
+    on(textbox.stored_string) do s
+        F[] = parse(Float64, s)
+    end
     # Update heatmap on slider movement
     on(slider.value) do i
         u = reshape(sol[i][1:N^2], N, N)
         heat_obs[] = u
     end
 
+
     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
+    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
+
+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[])")
+    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
+
+function build_ui(U, V, Du, Dv, F, k, params_obs, heat_obs)
+    reset!(U, V, heat_obs)
+    fig = Figure(size=(800, 800))
+    gh = GridLayout(fig[1, 1])
+    ax = Axis(gh[1, 1])
+    hm = Makie.heatmap!(ax, heat_obs, colormap=:viridis)
+    deactivate_interaction!(ax, :rectanglezoom)
+    ax.aspect = DataAspect()
+
+    run_label = Observable{String}("Run")
+    stepsize = Observable(30)
+    spoint = select_point(ax.scene)
+
+    # # 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
+
+
+    gh[1, 2] = textboxgrid = GridLayout(ax.scene, tellwidth=false)
+    rowsize!(gh, 1, Relative(1.0))
+
+
+    param_box!(textboxgrid, 1, "Du", Du)
+    param_box!(textboxgrid, 2, "Dv", Dv)
+    param_box!(textboxgrid, 3, "Feed", F)
+    param_box!(textboxgrid, 4, "kill", k)
+
+    # Timer and state for animation
+    running = Observable(false)
+
+    on(running) do r
+        run_label[] = r ? "Pause" : "Run"
+    end
+    on(speed_slider) do s
+        try
+            stepsize[] = s
+            println("Changed stepsize to $s")
+        catch
+            @warn "Invalid input for $s"
+        end
+    end
+    # Button Listeners
+    on(btn_step.clicks) do _
+        multi_step!((U, V), stepsize[], heat_obs, params_obs)
+    end
+
+    on(btn_start.clicks) do _
+        running[] = !running[]
+    end
+
+    on(btn_start.clicks) do _
+        @async while running[]
+            multi_step!((U, V), stepsize[], heat_obs, params_obs)
+            sleep(0.0015)  # ~20 FPS
+        end
+    end
+
+    on(btn_reset.clicks) do _
+        running[] = false
+        reset!(U, V, heat_obs)
+
+    end
+    on(spoint) do pt
+        r = 5
+        if pt === nothing
+            return
+        end
+        x, y = pt
+        println("params_obs[].N, ", params_obs[].N)
+        i, j = coord_to_index(x, y, params_obs[].N)
+
+        # get corners of square that will get filled with concentration
+        imin = max(i - r, 1)
+        imax = min(i + r, params_obs[].N)
+        jmin = max(j - r, 1)
+        jmax = min(j + r, params_obs[].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
+    return fig
+end
+
+
+export step_through_solution, build_ui
 end