diff --git a/scripts/main.jl b/scripts/main.jl
new file mode 100644
index 0000000..e5566c6
--- /dev/null
+++ b/scripts/main.jl
@@ -0,0 +1,28 @@
+include("../src/gray_scott_solver.jl")
+include("../src/visualization.jl")
+include("../src/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
index 82a719c..4748261 100644
--- a/scripts/run_simulation.jl
+++ b/scripts/run_simulation.jl
@@ -1,8 +1,10 @@
 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
@@ -13,12 +15,13 @@ F, k = Observable(0.0617), 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 current_params()
-    println(params_obs)
-end
 
-function update_params!(params_obs, u, v, feed, kill)
+
+function update_params!(params_obs::Observable, u, v, feed, kill)
     old = params_obs[]
     params_obs[] = GSParams(old.N, old.dx, u, v, feed, kill)
     current_params()
@@ -71,6 +74,7 @@ end
 # Observable holding current U for heatmap
 heat_obs = Observable(U)
 
+
 function laplacian5(f)
     h2 = dx^2
     left = f[2:end-1, 1:end-2]
@@ -132,14 +136,51 @@ ax.aspect = DataAspect()
 
 spoint = select_point(ax.scene)
 
-# # Controls
-run_label = Observable("Run")
-fig[2, 1] = buttongrid = GridLayout(tellwidth=false)
-btn_step = Button(buttongrid[1, 1], label="Step")
-btn_start = Button(buttongrid[1, 2], width=50, label=run_label)
-btn_reset = Button(buttongrid[1, 3], label="Reset")
+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
 
-gh[1, 2] = textboxgrid = GridLayout(tellwidth=false)
+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)
@@ -210,9 +251,10 @@ on(running) do r
     run_label[] = r ? "Pause" : "Run"
 end
 
+
 # Button Listeners
 on(btn_step.clicks) do _
-    multi_step!((U, V), 30)
+    multi_step!((U, V), stepsize[])
 end
 
 on(btn_start.clicks) do _
@@ -221,7 +263,7 @@ end
 
 on(btn_start.clicks) do _
     @async while running[]
-        multi_step!((U, V), 30)
+        multi_step!((U, V), stepsize[])
         sleep(0.05)  # ~20 FPS
     end
 end
diff --git a/src/gray_scott_solver.jl b/src/gray_scott_solver.jl
new file mode 100644
index 0000000..ca1840b
--- /dev/null
+++ b/src/gray_scott_solver.jl
@@ -0,0 +1,55 @@
+module GrayScottSolver
+using Observables
+include("constants.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
+
+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
index 259b5d8..3910c8f 100644
--- a/src/solver.jl
+++ b/src/solver.jl
@@ -1,5 +1,6 @@
 using DifferentialEquations
 using Random
+include("constants.jl")
 using .Constants
 
 """
diff --git a/src/visualization.jl b/src/visualization.jl
index 7f1514a..ec4e778 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,120 @@ 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)
+    # # 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
+    return fig
+end
+
+
+export step_through_solution, build_ui
 end