diff --git a/Manifest.toml b/Manifest.toml
index 23b9f57..6aaa7a8 100644
--- a/Manifest.toml
+++ b/Manifest.toml
@@ -2,7 +2,7 @@
 
 julia_version = "1.11.5"
 manifest_format = "2.0"
-project_hash = "0d147eaf78630b05e95b5317275a880443440272"
+project_hash = "0465e7f271b455d66e099c4ff05fca7f6e5b2b64"
 
 [[deps.ADTypes]]
 git-tree-sha1 = "e2478490447631aedba0823d4d7a80b2cc8cdb32"
diff --git a/Project.toml b/Project.toml
index 19c0af6..de40aee 100644
--- a/Project.toml
+++ b/Project.toml
@@ -1,4 +1,5 @@
 [deps]
 DifferentialEquations = "0c46a032-eb83-5123-abaf-570d42b7fbaa"
 GLMakie = "e9467ef8-e4e7-5192-8a1a-b1aee30e663a"
+Makie = "ee78f7c6-11fb-53f2-987a-cfe4a2b5a57a"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
diff --git a/scripts/run_simulation.jl b/scripts/run_simulation.jl
index 3e43bb4..6146ea8 100644
--- a/scripts/run_simulation.jl
+++ b/scripts/run_simulation.jl
@@ -4,8 +4,8 @@ using .AnimalFurFHN
 include("../src/visualization.jl")
 using .Visualization
 
-N = 100
-tspan = (0.0, 1000.0)
+N = 256
+tspan = (0.0, 5000.0)
 
 sol = AnimalFurFHN.run_simulation(tspan, N)
 
diff --git a/src/laplacian.jl b/src/laplacian.jl
index f9f4965..6a04d64 100644
--- a/src/laplacian.jl
+++ b/src/laplacian.jl
@@ -13,8 +13,9 @@
     - `Vector{Float64}`: A flattened (vectorized) representation of the approximated Laplacian values for each element in `U`. The boundary conditions are handled circularly.
 """
 function laplacian(U::Matrix{Float64}, N::Int, h::Float64)
-    # shifts matrices and sums them up
-    padded = circshift(U, (-1, 0)) .+ circshift(U, (1, 0)) .+
-             circshift(U, (0, -1)) .+ circshift(U, (0, 1)) .- 4 .* U
-    return vec(padded) ./ h^2
-end
+    Δ = zeros(N, N)
+    for i in 2:N-1, j in 2:N-1
+        Δ[i, j] = (U[i+1, j] + U[i-1, j] + U[i, j+1] + U[i, j-1] - 4 * U[i, j]) / h^2
+    end
+    return vec(Δ)
+end
\ No newline at end of file
diff --git a/src/solver.jl b/src/solver.jl
index 8999319..3ddd9f0 100644
--- a/src/solver.jl
+++ b/src/solver.jl
@@ -17,7 +17,7 @@ using .Constants
     # Returns
     - `du`: calculated derivatives put back into the du array
 """
-function fhn!(du, u, p::FHNParams, t = 0)
+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
 
@@ -30,6 +30,70 @@ 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.534522
+    end
+
+    safe_block!(u, p, p)
+
+    return vec(u), vec(v)
+end
+
+function column_ic(N)
+    u = fill(0.534522, N, N)
+    v = fill(0.381802, 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.534522
+
+    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.534522, N, N)
+    v = fill(0.381802, N, N)
+    cx, cy = div(N, 2), div(N, 2)  # center of matrix
+    radius = 0.250 * 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.534522
+        end
+    end
+
+    return vec(u), vec(v)
+end
+
+
 """
     run_simulation(tspan::Tuple{Float64,Float64}, N::Int)
 
@@ -44,18 +108,21 @@ end
 """
 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)
+    p = FHNParams(N=N, dx=1.0, Du=0.00016, Dv=0.001, ϵ=0.1, a=0.5, b=0.9)
 
     # 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))
+    #u0 = vec(0.4 .+ 0.01 .* rand(N, N))
+    #v0 = vec(0.4 .+ 0.01 .* rand(N, N))
 
-    y0 = vcat(u0, v0)
+    # Or use this
+    u0, v0 = center_band_ic(N)
+
+    y0 = vcat(vcat(u0, v0))
 
     prob = ODEProblem(fhn!, y0, tspan, p)
-    sol = solve(prob, BS3(), saveat=1.0)  # You can try `Rosenbrock23()` too
+    sol = solve(prob, BS3(), saveat=3.0)  # You can try `Rosenbrock23()` too
 
     return sol
 end
diff --git a/src/visualization.jl b/src/visualization.jl
index 7f1514a..00a0b71 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=:magma)
+    hmap = heatmap!(ax, heat_obs, colormap=:berlin)
 
     # Update heatmap on slider movement
     on(slider.value) do i