diff --git a/scripts/run_simulation.jl b/scripts/run_simulation.jl
index 0e253a3..b724908 100644
--- a/scripts/run_simulation.jl
+++ b/scripts/run_simulation.jl
@@ -5,7 +5,7 @@ include("../src/visualization.jl")
 using .Visualization
 
 N = 256
-tspan = (0.0, 1500.0)
+tspan = (1.0, 3000.0)
 
 sol = AnimalFurFHN.run_simulation(tspan, N)
 
diff --git a/src/AnimalFurFHN.jl b/src/AnimalFurFHN.jl
index 060abd7..f3f77c3 100644
--- a/src/AnimalFurFHN.jl
+++ b/src/AnimalFurFHN.jl
@@ -3,6 +3,7 @@ module AnimalFurFHN
 
 include("constants.jl")
 include("laplacian.jl")
+include("utils.jl")
 include("solver.jl")
 
 export run_simulation   # Make sure this is here!
diff --git a/src/solver.jl b/src/solver.jl
index a2be899..07df5c5 100644
--- a/src/solver.jl
+++ b/src/solver.jl
@@ -30,162 +30,6 @@ function fhn!(du, u, p::FHNParams, t=0)
     du .= vcat(vec(fu), vec(fv))
 end
 
-# helper functions for filling cells in specific places of the matrix
-function blocks_ic(N)
-    u = fill(1.0, N, N)
-    v = fill(0.0, N, N)
-    p = div(N, 2)
-
-    function safe_block!(u, row_center, col_center)
-        row_start = max(row_center - 8, 1)
-        row_end   = min(row_center + 7, N)
-        col_start = max(col_center - 8, 1)
-        col_end   = min(col_center + 7, N)
-        u[row_start:row_end, col_start:col_end] .= -0.01
-    end
-
-    safe_block!(u, p, p)
-
-    return vec(u), vec(v)
-end
-
-function column_ic(N)
-    u = fill(0.01, N, N)
-    v = fill(0.99, N, N)
-
-    col_center = div(N, 2)
-    col_width = 8  # You can adjust this
-
-    col_start = max(col_center - div(col_width, 2), 1)
-    col_end = min(col_center + div(col_width, 2) - 1, N)
-
-    u[col_start:col_end, :] .= -0.01
-
-    return vec(u), vec(v)
-end
-
-function two_rows_edge_distance_ic(N)
-    row_width = 8
-    distance_from_edge = 50
-    u = fill(0.01, N, N)
-    v = fill(0.99, N, N)
-
-    # --- Input Validation ---
-    if row_width <= 0 || distance_from_edge < 0
-        error("row_width must be positive and distance_from_edge must be non-negative.")
-    end
-
-    # Calculate column 1 (from the left edge)
-    col1_start = distance_from_edge + 1
-    col1_end = col1_start + row_width - 1
-
-    # Calculate column 2 (from the right edge)
-    col2_end = N - distance_from_edge
-    col2_start = col2_end - row_width + 1
-
-    # --- Further Validation for placement ---
-    if col1_end > N || col2_start < 1
-        error("Columns go out of bounds. Adjust N, row_width, or distance_from_edge.")
-    end
-    if col1_end >= col2_start # Check for overlap or touching
-        error("Columns overlap or touch. Adjust N, row_width, or distance_from_edge such that 2 * (distance_from_edge + row_width) <= N.")
-    end
-
-    # Apply the first column
-    u[:, col1_start:col1_end] .= -0.01
-
-    # Apply the second column
-    u[:, col2_start:col2_end] .= -0.01
-
-    return vec(u), vec(v)
-end
-
-function center_band_ic(N)
-    u = fill(0.0, N, N)
-    v = fill(0.0, N, N)
-
-    band_width = div(N, 8)
-    row_start = div(N, 2) - div(band_width, 2)
-    row_end = div(N, 2) + div(band_width, 2)
-
-    u[row_start:row_end, :] .= 0.1 .+ 0.01 .* randn(band_width + 1, N)
-    v[row_start:row_end, :] .= 0.1 .+ 0.01 .* randn(band_width + 1, N)
-
-    return vec(u), vec(v)
-end
-
-function circle_ic(N)
-    u = fill(0.01, N, N)
-    v = fill(0.99, N, N)
-    cx, cy = div(N, 2), div(N, 2)  # center of matrix
-    radius = 0.125 * N  # circle radius = 3/4 of N divided by 2
-
-    for i in 1:N, j in 1:N
-        if (i - cx)^2 + (j - cy)^2 ≤ radius^2
-            u[i, j] = -0.01
-        end
-    end
-
-    return vec(u), vec(v)
-end
-
-function three_circles_random_ic(N)
-    u = fill(0.01, N, N)
-    v = fill(0.99, N, N)
-    radius = 0.125 * N
-
-    # Define the bounds for random centers to ensure the circle stays within the matrix
-    min_coord = ceil(Int, radius) + 1
-    max_coord = floor(Int, N - radius)
-
-    if min_coord > max_coord
-        error("Matrix size N is too small to place circles of this radius without overlap or going out of bounds.")
-    end
-
-    for _ in 1:5 # Place 3 circles
-        # Generate random center coordinates
-        cx = rand(min_coord:max_coord)
-        cy = rand(min_coord:max_coord)
-
-        # Apply the circle to the matrix
-        for i in 1:N, j in 1:N
-            if (i - cx)^2 + (j - cy)^2 ≤ radius^2
-                u[i, j] = -0.01
-            end
-        end
-    end
-
-    return vec(u), vec(v)
-end
-
-function squiggle_ic(N, Lx=400.0, Ly=400.0)
-    uplus = 0.01
-    vplus = 0.99
-    uminus = -uplus
-
-    # Create coordinate grids
-    x = LinRange(0, Lx, N)
-    y = LinRange(0, Ly, N)
-    X = repeat(x', N, 1)  # Transposed to align with meshgrid X
-    Y = repeat(y, 1, N)   # Broadcasted to align with meshgrid Y
-
-    # Squiggle pattern
-    cos_term = 0.05 * Lx .* sin.(10 * 2π .* Y ./ Ly .+ π * 0.3)
-    exp_term = exp.(-((Y .- Ly / 2) ./ (0.1 * Ly)) .^ 2)
-    width = 0.05 * Lx
-    Z = exp.(-((X .- Lx / 2 .+ cos_term .* exp_term) ./ width) .^ 2)
-
-
-    u = fill(uplus, N, N)
-    v = fill(vplus, N, N)
-
-    # Apply squiggle
-    u[Z .> 0.8] .= uminus
-
-    return vec(u), vec(v)
-end
-
-
 """
     run_simulation(tspan::Tuple{Float64,Float64}, N::Int)
 
@@ -199,21 +43,11 @@ end
     - `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=0.016, Dv=0.1, ϵ=0.1, a=0.5, b=0.9)
+    # Initial conditions
+    # params, y0 = zebra_conditions(N) # tspan of ~3500 enough
+    params, y0 = cheetah_conditions(N)
 
-    # Initial conditions (random noise)
-    #Random.seed!(4321)
-    # Create two vectors with length N*N with numbers between 0.1 and 0.11
-    #u0 = vec(0.4 .+ 0.01 .* rand(N, N))
-    #v0 = vec(0.4 .+ 0.01 .* rand(N, N))
-
-    # Or use this
-    u0, v0 = two_rows_edge_distance_ic(N)
-
-    y0 = vcat(vcat(u0, v0))
-
-    prob = ODEProblem(fhn!, y0, tspan, p)
+    prob = ODEProblem(fhn!, y0, tspan, params)
     sol = solve(prob, BS3(), saveat=10.0)  # You can try `Rosenbrock23()` too
 
     return sol
diff --git a/src/utils.jl b/src/utils.jl
new file mode 100644
index 0000000..cf54588
--- /dev/null
+++ b/src/utils.jl
@@ -0,0 +1,184 @@
+
+# initial conditions for different patterns
+function zebra_conditions(N)
+    # Turing-spot parameters
+    params = FHNParams(N=N, dx=1.0, Du=0.016, Dv=0.1, ϵ=0.1, a=0.5, b=0.9)
+
+    # Or use this
+    u0, v0 = two_rows_edge_distance_ic(N)
+
+    return params, vcat(u0, v0)
+end
+
+function cheetah_conditions(N)
+    # Turing-spot parameters
+    params = FHNParams(N=N, dx=1.0, Du=0.0025, Dv=0.6, ϵ=0.05, a=0.7, b=0.8)
+
+    # Approximate Homogenous Steady State (HSS) for a=0.7, b=0.8
+    u_hss = -1.2
+    v_hss = -0.625
+
+    # Small perturbations around the HSS
+    # rand(N,N) generates values between 0 and 1.
+    # (2 .* rand(N,N) .- 1) generates values between -1 and 1 symmetrically around 0.
+    perturbation_amplitude = 0.01
+    u0 = vec(u_hss .+ perturbation_amplitude .* (2 .* rand(N, N) .- 1))
+    v0 = vec(v_hss .+ perturbation_amplitude .* (2 .* rand(N, N) .- 1))
+
+    return params, vcat(u0, v0)
+end
+
+# helper functions for filling cells in specific places of the matrix
+function blocks_ic(N)
+    u = fill(1.0, N, N)
+    v = fill(0.0, N, N)
+    p = div(N, 2)
+
+    function safe_block!(u, row_center, col_center)
+        row_start = max(row_center - 8, 1)
+        row_end   = min(row_center + 7, N)
+        col_start = max(col_center - 8, 1)
+        col_end   = min(col_center + 7, N)
+        u[row_start:row_end, col_start:col_end] .= -0.01
+    end
+
+    safe_block!(u, p, p)
+
+    return vec(u), vec(v)
+end
+
+function column_ic(N)
+    u = fill(0.01, N, N)
+    v = fill(0.99, N, N)
+
+    col_center = div(N, 2)
+    col_width = 8  # You can adjust this
+
+    col_start = max(col_center - div(col_width, 2), 1)
+    col_end = min(col_center + div(col_width, 2) - 1, N)
+
+    u[col_start:col_end, :] .= -0.01
+
+    return vec(u), vec(v)
+end
+
+function two_rows_edge_distance_ic(N)
+    row_width = 8
+    distance_from_edge = 50
+    u = fill(0.01, N, N)
+    v = fill(0.99, N, N)
+
+    # --- Input Validation ---
+    if row_width <= 0 || distance_from_edge < 0
+        error("row_width must be positive and distance_from_edge must be non-negative.")
+    end
+
+    # Calculate column 1 (from the left edge)
+    col1_start = distance_from_edge + 1
+    col1_end = col1_start + row_width - 1
+
+    # Calculate column 2 (from the right edge)
+    col2_end = N - distance_from_edge
+    col2_start = col2_end - row_width + 1
+
+    # --- Further Validation for placement ---
+    if col1_end > N || col2_start < 1
+        error("Columns go out of bounds. Adjust N, row_width, or distance_from_edge.")
+    end
+    if col1_end >= col2_start # Check for overlap or touching
+        error("Columns overlap or touch. Adjust N, row_width, or distance_from_edge such that 2 * (distance_from_edge + row_width) <= N.")
+    end
+
+    # Apply the first column
+    u[:, col1_start:col1_end] .= -0.01
+
+    # Apply the second column
+    u[:, col2_start:col2_end] .= -0.01
+
+    return vec(u), vec(v)
+end
+
+function center_band_ic(N)
+    u = fill(0.0, N, N)
+    v = fill(0.0, N, N)
+
+    band_width = div(N, 8)
+    row_start = div(N, 2) - div(band_width, 2)
+    row_end = div(N, 2) + div(band_width, 2)
+
+    u[row_start:row_end, :] .= 0.1 .+ 0.01 .* randn(band_width + 1, N)
+    v[row_start:row_end, :] .= 0.1 .+ 0.01 .* randn(band_width + 1, N)
+
+    return vec(u), vec(v)
+end
+
+function circle_ic(N)
+    u = fill(0.01, N, N)
+    v = fill(0.99, N, N)
+    cx, cy = div(N, 2), div(N, 2)  # center of matrix
+    radius = 0.125 * N  # circle radius = 3/4 of N divided by 2
+
+    for i in 1:N, j in 1:N
+        if (i - cx)^2 + (j - cy)^2 ≤ radius^2
+            u[i, j] = -0.01
+        end
+    end
+
+    return vec(u), vec(v)
+end
+
+function three_circles_random_ic(N)
+    u = fill(0.01, N, N)
+    v = fill(0.99, N, N)
+    radius = 0.125 * N
+
+    # Define the bounds for random centers to ensure the circle stays within the matrix
+    min_coord = ceil(Int, radius) + 1
+    max_coord = floor(Int, N - radius)
+
+    if min_coord > max_coord
+        error("Matrix size N is too small to place circles of this radius without overlap or going out of bounds.")
+    end
+
+    for _ in 1:5 # Place 3 circles
+        # Generate random center coordinates
+        cx = rand(min_coord:max_coord)
+        cy = rand(min_coord:max_coord)
+
+        # Apply the circle to the matrix
+        for i in 1:N, j in 1:N
+            if (i - cx)^2 + (j - cy)^2 ≤ radius^2
+                u[i, j] = -0.01
+            end
+        end
+    end
+
+    return vec(u), vec(v)
+end
+
+function squiggle_ic(N, Lx=400.0, Ly=400.0)
+    uplus = 0.01
+    vplus = 0.99
+    uminus = -uplus
+
+    # Create coordinate grids
+    x = LinRange(0, Lx, N)
+    y = LinRange(0, Ly, N)
+    X = repeat(x', N, 1)  # Transposed to align with meshgrid X
+    Y = repeat(y, 1, N)   # Broadcasted to align with meshgrid Y
+
+    # Squiggle pattern
+    cos_term = 0.05 * Lx .* sin.(10 * 2π .* Y ./ Ly .+ π * 0.3)
+    exp_term = exp.(-((Y .- Ly / 2) ./ (0.1 * Ly)) .^ 2)
+    width = 0.05 * Lx
+    Z = exp.(-((X .- Lx / 2 .+ cos_term .* exp_term) ./ width) .^ 2)
+
+
+    u = fill(uplus, N, N)
+    v = fill(vplus, N, N)
+
+    # Apply squiggle
+    u[Z .> 0.8] .= uminus
+
+    return vec(u), vec(v)
+end
diff --git a/src/visualization.jl b/src/visualization.jl
index 00a0b71..e73d9f1 100644
--- a/src/visualization.jl
+++ b/src/visualization.jl
@@ -21,7 +21,7 @@ function step_through_solution(sol, N::Int)
     # 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=:berlin)
+    hmap = heatmap!(ax, heat_obs, colormap=:RdGy)
 
     # Update heatmap on slider movement
     on(slider.value) do i