diff --git a/Manifest.toml b/Manifest.toml
index eceefdd..6c62bc0 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.4"
+julia_version = "1.11.5"
 manifest_format = "2.0"
-project_hash = "b5f5f0b50b1e0b7dd05bb93e0b1bb03fea235d53"
+project_hash = "5fd84347cd356de5b819aa3ea793fa63c45180e4"
 
 [[deps.ADTypes]]
 git-tree-sha1 = "e2478490447631aedba0823d4d7a80b2cc8cdb32"
@@ -194,9 +194,9 @@ version = "0.1.6"
 
 [[deps.BoundaryValueDiffEq]]
 deps = ["ADTypes", "BoundaryValueDiffEqAscher", "BoundaryValueDiffEqCore", "BoundaryValueDiffEqFIRK", "BoundaryValueDiffEqMIRK", "BoundaryValueDiffEqMIRKN", "BoundaryValueDiffEqShooting", "DiffEqBase", "FastClosures", "ForwardDiff", "LinearAlgebra", "Reexport", "SciMLBase"]
-git-tree-sha1 = "6606f3f7d43b038a8c34ee4cdc31bbd93b447407"
+git-tree-sha1 = "d6ec33e4516b2e790a64128afdb54f3b536667a7"
 uuid = "764a87c0-6b3e-53db-9096-fe964310641d"
-version = "5.17.0"
+version = "5.18.0"
 
     [deps.BoundaryValueDiffEq.extensions]
     BoundaryValueDiffEqODEInterfaceExt = "ODEInterface"
@@ -206,39 +206,39 @@ version = "5.17.0"
 
 [[deps.BoundaryValueDiffEqAscher]]
 deps = ["ADTypes", "AlmostBlockDiagonals", "BoundaryValueDiffEqCore", "ConcreteStructs", "DiffEqBase", "DifferentiationInterface", "FastClosures", "ForwardDiff", "LinearAlgebra", "PreallocationTools", "RecursiveArrayTools", "Reexport", "SciMLBase", "Setfield"]
-git-tree-sha1 = "add7560694794f30eec7803b9c4036efb96bf493"
+git-tree-sha1 = "64a777e06d995f677c86c7ddbb85f393074a0877"
 uuid = "7227322d-7511-4e07-9247-ad6ff830280e"
-version = "1.6.0"
+version = "1.7.0"
 
 [[deps.BoundaryValueDiffEqCore]]
 deps = ["ADTypes", "Adapt", "ArrayInterface", "ConcreteStructs", "DiffEqBase", "ForwardDiff", "LineSearch", "LinearAlgebra", "Logging", "NonlinearSolveFirstOrder", "PreallocationTools", "RecursiveArrayTools", "Reexport", "SciMLBase", "Setfield", "SparseArrays", "SparseConnectivityTracer", "SparseMatrixColorings"]
-git-tree-sha1 = "1e97d81b6ea5e8e938c57a3c640a5cff1e181a2f"
+git-tree-sha1 = "c83bf97da90dd379b1e3f4d9c6f3d0ae48eb0b29"
 uuid = "56b672f2-a5fe-4263-ab2d-da677488eb3a"
-version = "1.9.0"
+version = "1.10.0"
 
 [[deps.BoundaryValueDiffEqFIRK]]
 deps = ["ADTypes", "ArrayInterface", "BandedMatrices", "BoundaryValueDiffEqCore", "ConcreteStructs", "DiffEqBase", "DifferentiationInterface", "FastAlmostBandedMatrices", "FastClosures", "ForwardDiff", "LinearAlgebra", "PreallocationTools", "PrecompileTools", "Preferences", "RecursiveArrayTools", "Reexport", "SciMLBase", "Setfield", "SparseArrays"]
-git-tree-sha1 = "525c65da4b46271b4b5fc2178ccde16c99c21a41"
+git-tree-sha1 = "e58ee9acfc6dce6dcc368fc0483952bec5625513"
 uuid = "85d9eb09-370e-4000-bb32-543851f73618"
-version = "1.7.0"
+version = "1.8.1"
 
 [[deps.BoundaryValueDiffEqMIRK]]
 deps = ["ADTypes", "ArrayInterface", "BandedMatrices", "BoundaryValueDiffEqCore", "ConcreteStructs", "DiffEqBase", "DifferentiationInterface", "FastAlmostBandedMatrices", "FastClosures", "ForwardDiff", "LinearAlgebra", "PreallocationTools", "PrecompileTools", "Preferences", "RecursiveArrayTools", "Reexport", "SciMLBase", "Setfield", "SparseArrays"]
-git-tree-sha1 = "f729bdedaedb537bf7883c9f71e7577e1c7a07f6"
+git-tree-sha1 = "43debeee94167e2dc744f4a385213c4f0d16b4c3"
 uuid = "1a22d4ce-7765-49ea-b6f2-13c8438986a6"
-version = "1.7.0"
+version = "1.8.1"
 
 [[deps.BoundaryValueDiffEqMIRKN]]
 deps = ["ADTypes", "ArrayInterface", "BandedMatrices", "BoundaryValueDiffEqCore", "ConcreteStructs", "DiffEqBase", "DifferentiationInterface", "FastAlmostBandedMatrices", "FastClosures", "ForwardDiff", "LinearAlgebra", "PreallocationTools", "PrecompileTools", "Preferences", "RecursiveArrayTools", "Reexport", "SciMLBase", "Setfield", "SparseArrays"]
-git-tree-sha1 = "dfde7afe34137182a4bb8cd98f01d870772be999"
+git-tree-sha1 = "1d92c9f7567b627514e143a3caf93af6d235c2db"
 uuid = "9255f1d6-53bf-473e-b6bd-23f1ff009da4"
-version = "1.6.0"
+version = "1.7.1"
 
 [[deps.BoundaryValueDiffEqShooting]]
 deps = ["ADTypes", "ArrayInterface", "BandedMatrices", "BoundaryValueDiffEqCore", "ConcreteStructs", "DiffEqBase", "DifferentiationInterface", "FastAlmostBandedMatrices", "FastClosures", "ForwardDiff", "LinearAlgebra", "PreallocationTools", "PrecompileTools", "Preferences", "RecursiveArrayTools", "Reexport", "SciMLBase", "Setfield", "SparseArrays"]
-git-tree-sha1 = "255c15c2d262ea8d3c5db1e0f130f829362c0fd4"
+git-tree-sha1 = "1460c449c91c9bc4c8fb08fb5e878811ab38d667"
 uuid = "ed55bfe0-3725-4db6-871e-a1dc9f42a757"
-version = "1.7.0"
+version = "1.8.1"
 
 [[deps.BracketingNonlinearSolve]]
 deps = ["CommonSolve", "ConcreteStructs", "NonlinearSolveBase", "PrecompileTools", "Reexport", "SciMLBase"]
@@ -390,9 +390,9 @@ uuid = "2569d6c7-a4a2-43d3-a901-331e8e4be471"
 version = "0.2.3"
 
 [[deps.ConstructionBase]]
-git-tree-sha1 = "76219f1ed5771adbb096743bff43fb5fdd4c1157"
+git-tree-sha1 = "b4b092499347b18a015186eae3042f72267106cb"
 uuid = "187b0558-2788-49d3-abe0-74a17ed4e7c9"
-version = "1.5.8"
+version = "1.6.0"
 weakdeps = ["IntervalSets", "LinearAlgebra", "StaticArrays"]
 
     [deps.ConstructionBase.extensions]
@@ -457,9 +457,9 @@ version = "5.53.1"
 
 [[deps.DiffEqBase]]
 deps = ["ArrayInterface", "ConcreteStructs", "DataStructures", "DocStringExtensions", "EnumX", "EnzymeCore", "FastBroadcast", "FastClosures", "FastPower", "FunctionWrappers", "FunctionWrappersWrappers", "LinearAlgebra", "Logging", "Markdown", "MuladdMacro", "Parameters", "PrecompileTools", "Printf", "RecursiveArrayTools", "Reexport", "SciMLBase", "SciMLOperators", "SciMLStructures", "Setfield", "Static", "StaticArraysCore", "Statistics", "SymbolicIndexingInterface", "TruncatedStacktraces"]
-git-tree-sha1 = "a0e5b5669df9465bc3dd32ea4a8ddeefbc0f7b5c"
+git-tree-sha1 = "2d87d7bd165c1ca0d11923a9fabe90a9d71e88a6"
 uuid = "2b5f629d-d688-5b77-993f-72d75c75574e"
-version = "6.175.0"
+version = "6.176.0"
 
     [deps.DiffEqBase.extensions]
     DiffEqBaseCUDAExt = "CUDA"
@@ -533,9 +533,9 @@ version = "7.16.1"
 
 [[deps.DifferentiationInterface]]
 deps = ["ADTypes", "LinearAlgebra"]
-git-tree-sha1 = "c8d85ecfcbaef899308706bebdd8b00107f3fb43"
+git-tree-sha1 = "210933c93f39f832d92f9efbbe69a49c453db36d"
 uuid = "a0c0ee7d-e4b9-4e03-894e-1c5f64a51d63"
-version = "0.6.54"
+version = "0.7.1"
 
     [deps.DifferentiationInterface.extensions]
     DifferentiationInterfaceChainRulesCoreExt = "ChainRulesCore"
@@ -614,9 +614,9 @@ version = "0.25.120"
     Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [[deps.DocStringExtensions]]
-git-tree-sha1 = "e7b7e6f178525d17c720ab9c081e4ef04429f860"
+git-tree-sha1 = "7442a5dfe1ebb773c29cc2962a8980f47221d76c"
 uuid = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae"
-version = "0.9.4"
+version = "0.9.5"
 
 [[deps.Downloads]]
 deps = ["ArgTools", "FileWatching", "LibCURL", "NetworkOptions"]
@@ -635,9 +635,9 @@ uuid = "4e289a0a-7415-4d19-859d-a7e5c4648b56"
 version = "1.0.5"
 
 [[deps.EnzymeCore]]
-git-tree-sha1 = "1eb59f40a772d0fbd4cb75e00b3fa7f5f79c975a"
+git-tree-sha1 = "7d7822a643c33bbff4eab9c87ca8459d7c688db0"
 uuid = "f151be2c-9106-41f4-ab19-57ee4f262869"
-version = "0.8.9"
+version = "0.8.11"
 weakdeps = ["Adapt"]
 
     [deps.EnzymeCore.extensions]
@@ -694,9 +694,9 @@ version = "6.1.2+0"
 
 [[deps.FFTW]]
 deps = ["AbstractFFTs", "FFTW_jll", "LinearAlgebra", "MKL_jll", "Preferences", "Reexport"]
-git-tree-sha1 = "7de7c78d681078f027389e067864a8d53bd7c3c9"
+git-tree-sha1 = "797762812ed063b9b94f6cc7742bc8883bb5e69e"
 uuid = "7a1cc6ca-52ef-59f5-83cd-3a7055c09341"
-version = "1.8.1"
+version = "1.9.0"
 
 [[deps.FFTW_jll]]
 deps = ["Artifacts", "JLLWrappers", "Libdl"]
@@ -852,9 +852,9 @@ version = "2.13.4+0"
 
 [[deps.FreeTypeAbstraction]]
 deps = ["BaseDirs", "ColorVectorSpace", "Colors", "FreeType", "GeometryBasics", "Mmap"]
-git-tree-sha1 = "eaca92bac73aa42f68c57d1b8df1b746eeb2bdaa"
+git-tree-sha1 = "4ebb930ef4a43817991ba35db6317a05e59abd11"
 uuid = "663a7486-cb36-511b-a19d-713bb74d65c9"
-version = "0.10.7"
+version = "0.10.8"
 
 [[deps.FriBidi_jll]]
 deps = ["Artifacts", "JLLWrappers", "Libdl"]
@@ -898,9 +898,9 @@ version = "3.4.0+2"
 
 [[deps.GLMakie]]
 deps = ["ColorTypes", "Colors", "FileIO", "FixedPointNumbers", "FreeTypeAbstraction", "GLFW", "GeometryBasics", "LinearAlgebra", "Makie", "Markdown", "MeshIO", "ModernGL", "Observables", "PrecompileTools", "Printf", "ShaderAbstractions", "StaticArrays"]
-git-tree-sha1 = "0650df73822ce8808dc473e3e1c7581bdb367083"
+git-tree-sha1 = "aacae8d91f7da30979e189e8046c5bd655421687"
 uuid = "e9467ef8-e4e7-5192-8a1a-b1aee30e663a"
-version = "0.11.8"
+version = "0.12.0"
 
 [[deps.GPUArraysCore]]
 deps = ["Adapt"]
@@ -956,10 +956,10 @@ uuid = "3b182d85-2403-5c21-9c21-1e1f0cc25472"
 version = "1.3.15+0"
 
 [[deps.Graphs]]
-deps = ["ArnoldiMethod", "Compat", "DataStructures", "Distributed", "Inflate", "LinearAlgebra", "Random", "SharedArrays", "SimpleTraits", "SparseArrays", "Statistics"]
-git-tree-sha1 = "3169fd3440a02f35e549728b0890904cfd4ae58a"
+deps = ["ArnoldiMethod", "DataStructures", "Distributed", "Inflate", "LinearAlgebra", "Random", "SharedArrays", "SimpleTraits", "SparseArrays", "Statistics"]
+git-tree-sha1 = "c5abfa0ae0aaee162a3fbb053c13ecda39be545b"
 uuid = "86223c79-3864-5bf0-83f7-82e725a168b6"
-version = "1.12.1"
+version = "1.13.0"
 
 [[deps.GridLayoutBase]]
 deps = ["GeometryBasics", "InteractiveUtils", "Observables"]
@@ -1145,9 +1145,9 @@ version = "3.1.1+0"
 
 [[deps.JumpProcesses]]
 deps = ["ArrayInterface", "DataStructures", "DiffEqBase", "DiffEqCallbacks", "DocStringExtensions", "FunctionWrappers", "Graphs", "LinearAlgebra", "Markdown", "PoissonRandom", "Random", "RandomNumbers", "RecursiveArrayTools", "Reexport", "SciMLBase", "Setfield", "StaticArrays", "SymbolicIndexingInterface", "UnPack"]
-git-tree-sha1 = "216c196df09c8b80a40a2befcb95760eb979bcfd"
+git-tree-sha1 = "fb7fd516de38db80f50fe15e57d44da2836365e7"
 uuid = "ccbc3e58-028d-4f4c-8cd5-9ae44345cda5"
-version = "9.15.0"
+version = "9.16.0"
 weakdeps = ["FastBroadcast"]
 
 [[deps.KernelDensity]]
@@ -1308,9 +1308,9 @@ weakdeps = ["LineSearches"]
 
 [[deps.LineSearches]]
 deps = ["LinearAlgebra", "NLSolversBase", "NaNMath", "Parameters", "Printf"]
-git-tree-sha1 = "e4c3be53733db1051cc15ecf573b1042b3a712a1"
+git-tree-sha1 = "4adee99b7262ad2a1a4bbbc59d993d24e55ea96f"
 uuid = "d3d80556-e9d4-5f37-9878-2ab0fcc64255"
-version = "7.3.0"
+version = "7.4.0"
 
 [[deps.LinearAlgebra]]
 deps = ["Libdl", "OpenBLAS_jll", "libblastrampoline_jll"]
@@ -1319,9 +1319,9 @@ version = "1.11.0"
 
 [[deps.LinearSolve]]
 deps = ["ArrayInterface", "ChainRulesCore", "ConcreteStructs", "DocStringExtensions", "EnumX", "GPUArraysCore", "InteractiveUtils", "Krylov", "LazyArrays", "Libdl", "LinearAlgebra", "MKL_jll", "Markdown", "PrecompileTools", "Preferences", "RecursiveArrayTools", "Reexport", "SciMLBase", "SciMLOperators", "Setfield", "StaticArraysCore", "UnPack"]
-git-tree-sha1 = "c618a6a774d5712c6bf02dbcceb51b6dc6b9bb89"
+git-tree-sha1 = "c0d1a91a50af6778863d320761f807f641f74935"
 uuid = "7ed4a6bd-45f5-4d41-b270-4a48e9bafcae"
-version = "3.16.0"
+version = "3.17.0"
 
     [deps.LinearSolve.extensions]
     LinearSolveBandedMatricesExt = "BandedMatrices"
@@ -1392,15 +1392,15 @@ version = "0.5.16"
 
 [[deps.Makie]]
 deps = ["Animations", "Base64", "CRC32c", "ColorBrewer", "ColorSchemes", "ColorTypes", "Colors", "Contour", "Dates", "DelaunayTriangulation", "Distributions", "DocStringExtensions", "Downloads", "FFMPEG_jll", "FileIO", "FilePaths", "FixedPointNumbers", "Format", "FreeType", "FreeTypeAbstraction", "GeometryBasics", "GridLayoutBase", "ImageBase", "ImageIO", "InteractiveUtils", "Interpolations", "IntervalSets", "InverseFunctions", "Isoband", "KernelDensity", "LaTeXStrings", "LinearAlgebra", "MacroTools", "MakieCore", "Markdown", "MathTeXEngine", "Observables", "OffsetArrays", "PNGFiles", "Packing", "PlotUtils", "PolygonOps", "PrecompileTools", "Printf", "REPL", "Random", "RelocatableFolders", "Scratch", "ShaderAbstractions", "Showoff", "SignedDistanceFields", "SparseArrays", "Statistics", "StatsBase", "StatsFuns", "StructArrays", "TriplotBase", "UnicodeFun", "Unitful"]
-git-tree-sha1 = "968f03dc65c8144a728f988051a88c78fe625e26"
+git-tree-sha1 = "29c5cf33287afa8c06711fec1b33d074f88dd354"
 uuid = "ee78f7c6-11fb-53f2-987a-cfe4a2b5a57a"
-version = "0.22.7"
+version = "0.23.0"
 
 [[deps.MakieCore]]
 deps = ["ColorTypes", "GeometryBasics", "IntervalSets", "Observables"]
-git-tree-sha1 = "733d910c70805e7114c82508bae99c6cdf004466"
+git-tree-sha1 = "6ec4aba99ce1a8c247822ce346b273d0d686e141"
 uuid = "20f20a25-4f0e-4fdf-b5d1-57303727442b"
-version = "0.9.3"
+version = "0.10.0"
 
 [[deps.ManualMemory]]
 git-tree-sha1 = "bcaef4fc7a0cfe2cba636d84cda54b5e4e4ca3cd"
@@ -1493,9 +1493,9 @@ version = "0.2.4"
 
 [[deps.NLSolversBase]]
 deps = ["ADTypes", "DifferentiationInterface", "Distributed", "FiniteDiff", "ForwardDiff"]
-git-tree-sha1 = "b14c7be6046e7d48e9063a0053f95ee0fc954176"
+git-tree-sha1 = "25a6638571a902ecfb1ae2a18fc1575f86b1d4df"
 uuid = "d41bc354-129a-5804-8e4c-c37616107c6c"
-version = "7.9.1"
+version = "7.10.0"
 
 [[deps.NLsolve]]
 deps = ["Distances", "LineSearches", "LinearAlgebra", "NLSolversBase", "Printf", "Reexport"]
@@ -1639,7 +1639,7 @@ version = "3.2.4+0"
 [[deps.OpenLibm_jll]]
 deps = ["Artifacts", "Libdl"]
 uuid = "05823500-19ac-5b8b-9628-191a04bc5112"
-version = "0.8.1+4"
+version = "0.8.5+0"
 
 [[deps.OpenSSL_jll]]
 deps = ["Artifacts", "JLLWrappers", "Libdl"]
@@ -1655,9 +1655,9 @@ version = "0.5.6+0"
 
 [[deps.Optim]]
 deps = ["Compat", "EnumX", "FillArrays", "ForwardDiff", "LineSearches", "LinearAlgebra", "NLSolversBase", "NaNMath", "PositiveFactorizations", "Printf", "SparseArrays", "StatsBase"]
-git-tree-sha1 = "31b3b1b8e83ef9f1d50d74f1dd5f19a37a304a1f"
+git-tree-sha1 = "61942645c38dd2b5b78e2082c9b51ab315315d10"
 uuid = "429524aa-4258-5aef-a3af-852621145aeb"
-version = "1.12.0"
+version = "1.13.2"
 
     [deps.Optim.extensions]
     OptimMOIExt = "MathOptInterface"
@@ -1690,15 +1690,15 @@ version = "1.2.0"
 
 [[deps.OrdinaryDiffEqBDF]]
 deps = ["ADTypes", "ArrayInterface", "DiffEqBase", "FastBroadcast", "LinearAlgebra", "MacroTools", "MuladdMacro", "OrdinaryDiffEqCore", "OrdinaryDiffEqDifferentiation", "OrdinaryDiffEqNonlinearSolve", "OrdinaryDiffEqSDIRK", "PrecompileTools", "Preferences", "RecursiveArrayTools", "Reexport", "StaticArrays", "TruncatedStacktraces"]
-git-tree-sha1 = "42755bd13fe56e9d9ce1bc005f8b206a6b56b731"
+git-tree-sha1 = "9124a686af119063bb4d3a8f87044a8f312fcad9"
 uuid = "6ad6398a-0878-4a85-9266-38940aa047c8"
-version = "1.5.1"
+version = "1.6.0"
 
 [[deps.OrdinaryDiffEqCore]]
 deps = ["ADTypes", "Accessors", "Adapt", "ArrayInterface", "DataStructures", "DiffEqBase", "DocStringExtensions", "EnumX", "FastBroadcast", "FastClosures", "FastPower", "FillArrays", "FunctionWrappersWrappers", "InteractiveUtils", "LinearAlgebra", "Logging", "MacroTools", "MuladdMacro", "Polyester", "PrecompileTools", "Preferences", "RecursiveArrayTools", "Reexport", "SciMLBase", "SciMLOperators", "SciMLStructures", "SimpleUnPack", "Static", "StaticArrayInterface", "StaticArraysCore", "SymbolicIndexingInterface", "TruncatedStacktraces"]
-git-tree-sha1 = "d29adfeb720dd7c251b216d91c4bd4fe67c087df"
+git-tree-sha1 = "08dac9c6672a4548439048089bac293759a897fd"
 uuid = "bbf590c4-e513-4bbe-9b18-05decba2e5d8"
-version = "1.26.0"
+version = "1.26.1"
 weakdeps = ["EnzymeCore"]
 
     [deps.OrdinaryDiffEqCore.extensions]
@@ -1712,9 +1712,9 @@ version = "1.4.0"
 
 [[deps.OrdinaryDiffEqDifferentiation]]
 deps = ["ADTypes", "ArrayInterface", "ConcreteStructs", "ConstructionBase", "DiffEqBase", "DifferentiationInterface", "FastBroadcast", "FiniteDiff", "ForwardDiff", "FunctionWrappersWrappers", "LinearAlgebra", "LinearSolve", "OrdinaryDiffEqCore", "SciMLBase", "SciMLOperators", "SparseArrays", "SparseMatrixColorings", "StaticArrayInterface", "StaticArrays"]
-git-tree-sha1 = "c78060115fa4ea9d70ac47fa49496acbc630aefa"
+git-tree-sha1 = "efecf0c4cc44e16251b0e718f08b0876b2a82b80"
 uuid = "4302a76b-040a-498a-8c04-15b101fed76b"
-version = "1.9.1"
+version = "1.10.0"
 
 [[deps.OrdinaryDiffEqExplicitRK]]
 deps = ["DiffEqBase", "FastBroadcast", "LinearAlgebra", "MuladdMacro", "OrdinaryDiffEqCore", "RecursiveArrayTools", "Reexport", "TruncatedStacktraces"]
@@ -1820,9 +1820,9 @@ version = "1.1.0"
 
 [[deps.OrdinaryDiffEqRosenbrock]]
 deps = ["ADTypes", "DiffEqBase", "DifferentiationInterface", "FastBroadcast", "FiniteDiff", "ForwardDiff", "LinearAlgebra", "LinearSolve", "MacroTools", "MuladdMacro", "OrdinaryDiffEqCore", "OrdinaryDiffEqDifferentiation", "Polyester", "PrecompileTools", "Preferences", "RecursiveArrayTools", "Reexport", "Static"]
-git-tree-sha1 = "063e5ff1447b3869856ed264b6dcbb21e6e8bdb0"
+git-tree-sha1 = "1ce0096d920e95773220e818f29bf4b37ea2bb78"
 uuid = "43230ef6-c299-4910-a778-202eb28ce4ce"
-version = "1.10.1"
+version = "1.11.0"
 
 [[deps.OrdinaryDiffEqSDIRK]]
 deps = ["ADTypes", "DiffEqBase", "FastBroadcast", "LinearAlgebra", "MacroTools", "MuladdMacro", "OrdinaryDiffEqCore", "OrdinaryDiffEqDifferentiation", "OrdinaryDiffEqNonlinearSolve", "RecursiveArrayTools", "Reexport", "SciMLBase", "TruncatedStacktraces"]
@@ -2169,9 +2169,9 @@ version = "0.1.0"
 
 [[deps.SciMLBase]]
 deps = ["ADTypes", "Accessors", "Adapt", "ArrayInterface", "CommonSolve", "ConstructionBase", "Distributed", "DocStringExtensions", "EnumX", "FunctionWrappersWrappers", "IteratorInterfaceExtensions", "LinearAlgebra", "Logging", "Markdown", "Moshi", "PrecompileTools", "Preferences", "Printf", "RecipesBase", "RecursiveArrayTools", "Reexport", "RuntimeGeneratedFunctions", "SciMLOperators", "SciMLStructures", "StaticArraysCore", "Statistics", "SymbolicIndexingInterface"]
-git-tree-sha1 = "ce947672206f6a3a2bee1017c690cfd5fd82d897"
+git-tree-sha1 = "04bbcdc8d1f7d6f667f75fbcc68728231e21fabe"
 uuid = "0bca4576-84f4-4d90-8ffe-ffa030f20462"
-version = "2.96.0"
+version = "2.101.0"
 
     [deps.SciMLBase.extensions]
     SciMLBaseChainRulesCoreExt = "ChainRulesCore"
@@ -2202,9 +2202,9 @@ version = "0.1.6"
 
 [[deps.SciMLOperators]]
 deps = ["Accessors", "ArrayInterface", "DocStringExtensions", "LinearAlgebra", "MacroTools"]
-git-tree-sha1 = "85608e4aaf758547ffc4030c908318b432114ec9"
+git-tree-sha1 = "3249fe77f322fe539e935ecb388c8290cd38a3fc"
 uuid = "c0aeaf25-5076-4817-a8d5-81caf7dfa961"
-version = "1.3.0"
+version = "1.3.1"
 weakdeps = ["SparseArrays", "StaticArraysCore"]
 
     [deps.SciMLOperators.extensions]
@@ -2308,9 +2308,9 @@ version = "1.11.0"
 
 [[deps.SparseConnectivityTracer]]
 deps = ["ADTypes", "DocStringExtensions", "FillArrays", "LinearAlgebra", "Random", "SparseArrays"]
-git-tree-sha1 = "fadb2d7010dd92912e5eb31a493613ad4b8c9583"
+git-tree-sha1 = "affde0bfd920cfcaa0944d3c0eb3a573fa9c4d1e"
 uuid = "9f842d2f-2579-4b1d-911e-f412cf18a3f5"
-version = "0.6.18"
+version = "0.6.20"
 
     [deps.SparseConnectivityTracer.extensions]
     SparseConnectivityTracerDataInterpolationsExt = "DataInterpolations"
@@ -2517,9 +2517,9 @@ version = "1.0.1"
 
 [[deps.Tables]]
 deps = ["DataAPI", "DataValueInterfaces", "IteratorInterfaceExtensions", "OrderedCollections", "TableTraits"]
-git-tree-sha1 = "598cd7c1f68d1e205689b1c2fe65a9f85846f297"
+git-tree-sha1 = "f2c1efbc8f3a609aadf318094f8fc5204bdaf344"
 uuid = "bd369af6-aec1-5ad0-b16a-f7cc5008161c"
-version = "1.12.0"
+version = "1.12.1"
 
 [[deps.Tar]]
 deps = ["ArgTools", "SHA"]
@@ -2539,15 +2539,15 @@ version = "1.11.0"
 
 [[deps.ThreadingUtilities]]
 deps = ["ManualMemory"]
-git-tree-sha1 = "2d529b6b22791f3e22e7ec5c60b9016e78f5f6bf"
+git-tree-sha1 = "d969183d3d244b6c33796b5ed01ab97328f2db85"
 uuid = "8290d209-cae3-49c0-8002-c8c24d57dab5"
-version = "0.5.4"
+version = "0.5.5"
 
 [[deps.TiffImages]]
-deps = ["ColorTypes", "DataStructures", "DocStringExtensions", "FileIO", "FixedPointNumbers", "IndirectArrays", "Inflate", "Mmap", "OffsetArrays", "PkgVersion", "ProgressMeter", "SIMD", "UUIDs"]
-git-tree-sha1 = "f21231b166166bebc73b99cea236071eb047525b"
+deps = ["ColorTypes", "DataStructures", "DocStringExtensions", "FileIO", "FixedPointNumbers", "IndirectArrays", "Inflate", "Mmap", "OffsetArrays", "PkgVersion", "PrecompileTools", "ProgressMeter", "SIMD", "UUIDs"]
+git-tree-sha1 = "02aca429c9885d1109e58f400c333521c13d48a0"
 uuid = "731e570b-9d59-4bfa-96dc-6df516fadf69"
-version = "0.11.3"
+version = "0.11.4"
 
 [[deps.TimerOutputs]]
 deps = ["ExprTools", "Printf"]
@@ -2599,14 +2599,16 @@ version = "0.4.1"
 
 [[deps.Unitful]]
 deps = ["Dates", "LinearAlgebra", "Random"]
-git-tree-sha1 = "d62610ec45e4efeabf7032d67de2ffdea8344bed"
+git-tree-sha1 = "d2282232f8a4d71f79e85dc4dd45e5b12a6297fb"
 uuid = "1986cc42-f94f-5a68-af5c-568840ba703d"
-version = "1.22.1"
-weakdeps = ["ConstructionBase", "InverseFunctions"]
+version = "1.23.1"
+weakdeps = ["ConstructionBase", "ForwardDiff", "InverseFunctions", "Printf"]
 
     [deps.Unitful.extensions]
     ConstructionBaseUnitfulExt = "ConstructionBase"
+    ForwardDiffExt = "ForwardDiff"
     InverseFunctionsUnitfulExt = "InverseFunctions"
+    PrintfExt = "Printf"
 
 [[deps.Wayland_jll]]
 deps = ["Artifacts", "EpollShim_jll", "Expat_jll", "JLLWrappers", "Libdl", "Libffi_jll", "XML2_jll"]
@@ -2615,10 +2617,10 @@ uuid = "a2964d1f-97da-50d4-b82a-358c7fce9d89"
 version = "1.23.1+0"
 
 [[deps.Wayland_protocols_jll]]
-deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
-git-tree-sha1 = "5db3e9d307d32baba7067b13fc7b5aa6edd4a19a"
+deps = ["Artifacts", "JLLWrappers", "Libdl"]
+git-tree-sha1 = "54b8a029ac145ebe8299463447fd1590b2b1d92f"
 uuid = "2381bf8a-dfd0-557d-9999-79630e7b1b91"
-version = "1.36.0+0"
+version = "1.44.0+0"
 
 [[deps.WebP]]
 deps = ["CEnum", "ColorTypes", "FileIO", "FixedPointNumbers", "ImageCore", "libwebp_jll"]
diff --git a/scripts/main.jl b/scripts/main.jl
index e38e373..9543a01 100644
--- a/scripts/main.jl
+++ b/scripts/main.jl
@@ -1,28 +1,58 @@
+include("../src/utils/constants.jl")
+#include("../src/fhn_solver.jl")
 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
+using .Constants
+using .Visualization
+
+# GSParams
+N = 128
+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[]))
+param_observables = (
+    Du=Du,
+    Dv=Dv,
+    F=F,
+    k=k,
+)
+
+"""
+# FHNParams
+N = 128
+dx = 1.0
+Du, Dv = Observable(0.016), Observable(0.1)
+ϵ, a, b = Observable(0.1), Observable(0.5), Observable(0.9)
+param_observables = (
+    Du=Du,
+    Dv=Dv,
+    ϵ=ϵ,
+    a=a,
+    b=b
+)
+"""
+params_obs = Observable(Constants.GSParams(N, dx, Du[], Dv[], F[], k[]))
+#params_obs = Observable(FHNParams(N=N, dx=dx, Du=Du[], Dv=Dv[], ϵ=ϵ[], a=a[], b=b[]))
+
+
 lift(Du, Dv, F, k) do u, v, f, ki
-    params_obs[] = GSParams(N, dx, u, v, f, ki)
+   params_obs[] = GSParams(N, dx, u, v, f, ki)
 end
 
+"""
+lift(Du, Dv, ϵ, a, b) do d_u, d_v, eps, aa, bb
+    params_obs[] = FHNParams(N=N, dx=dx, Du=d_u, Dv=d_v, ϵ=eps, a=aa, b=bb)
+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)
+fig = build_ui(U, V, param_observables, params_obs, heat_obs)
 display(fig)
diff --git a/scripts/show_frame.jl b/scripts/show_frame.jl
deleted file mode 100644
index e645396..0000000
--- a/scripts/show_frame.jl
+++ /dev/null
@@ -1,21 +0,0 @@
-using GLMakie
-include("../src/AnimalFurFHN.jl")
-using .AnimalFurFHN
-
-# === Parameters ===
-N = 100
-tspan = (0.0, 100.0)
-frame_index = 1  # final time step
-
-# === Run simulation ===
-sol = run_simulation(N, tspan)
-
-# === Extract activator u ===
-u = sol.u[frame_index][1:N^2]
-u_mat = reshape(u, N, N)
-
-# === Plot ===
-fig = Figure()
-ax = Axis(fig[1, 1])
-heatmap!(ax, u_mat, colormap=:viridis)
-fig
diff --git a/src/fhn_solver.jl b/src/fhn_solver.jl
new file mode 100644
index 0000000..75ca9b3
--- /dev/null
+++ b/src/fhn_solver.jl
@@ -0,0 +1,75 @@
+include("utils/laplacian.jl")
+
+using DifferentialEquations
+using Random
+using Observables
+
+using ..Constants
+using .Laplacian
+
+"""
+    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)
+    N = p.N
+    u_mat = reshape(u[1:N^2], N, N)
+    v_mat = reshape(u[N^2+1:end], N, N)
+
+    Δu = laplacian(u_mat, p.dx)
+    Δv = laplacian(v_mat, p.dx)
+
+    u_in = u_mat[2:end-1, 2:end-1]
+    v_in = v_mat[2:end-1, 2:end-1]
+
+    fu = p.Du * Δu .+ u_in .- u_in .^ 3 ./ 3 .- v_in
+    fv = p.Dv * Δv .+ p.ϵ * (u_in .+ p.a .- p.b .* v_in)
+
+    # Construct output with zero boundary (padding)
+    fu_full = zeros(N, N)
+    fv_full = zeros(N, N)
+    fu_full[2:end-1, 2:end-1] .= fu
+    fv_full[2:end-1, 2:end-1] .= fv
+
+    du .= vcat(vec(fu_full), vec(fv_full))
+end
+
+function step!(U, V, params_obs::Observable; dx=1)
+    p = params_obs[]
+
+    # Flatten initial condition (activation u, recovery v)
+    u0 = vec(U)
+    v0 = vec(V)
+    u_init = vcat(u0, v0)
+
+    # Define one integration step using your fhn! function
+    prob = ODEProblem((du, u, p, t) -> fhn!(du, u, p, t), u_init, (0.0, 0.1), p)
+    sol = solve(prob, Tsit5(); dt=0.1, save_start=false, saveat=10.0)
+
+    # Extract solution and reshape
+    u_new = reshape(sol[end][1:p.N^2], p.N, p.N)
+    v_new = reshape(sol[end][p.N^2+1:end], p.N, p.N)
+
+    # Update matrices in-place
+    U .= u_new
+    V .= v_new
+
+    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
diff --git a/src/gray_scott_solver.jl b/src/gray_scott_solver.jl
index 3668d26..e9e2c5f 100644
--- a/src/gray_scott_solver.jl
+++ b/src/gray_scott_solver.jl
@@ -1,30 +1,35 @@
-module GrayScottSolver
-using Observables
 include("utils/constants.jl")
 include("utils/laplacian.jl")
+
+using Observables
+
 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
+    # Extract parameters
+    p = params_obs[]
+    Du, Dv = p.Du, p.Dv
+    F, k = p.F, p.k
+
+    # Compute Laplacians on the interior
+    lap_u = laplacian(U, dx)  # shape (N-2, N-2)
+    lap_v = laplacian(V, dx)
+
+    # Extract interior values
     u = U[2:end-1, 2:end-1]
     v = V[2:end-1, 2:end-1]
 
+    # Gray-Scott update equations (Euler-style)
     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)
+    u_new = u .+ (Du .* lap_u .- uvv .+ F .* (1 .- u))
+    v_new = v .+ (Dv .* lap_v .+ uvv .- (F .+ k) .* v)
 
-    # Update with new values
+    # Update interior
     U[2:end-1, 2:end-1] .= u_new
     V[2:end-1, 2:end-1] .= v_new
 
-    # Periodic boundary conditions
+    # Apply periodic boundary conditions
     U[1, :] .= U[end-1, :]
     U[end, :] .= U[2, :]
     U[:, 1] .= U[:, end-1]
@@ -35,12 +40,10 @@ function step!(U, V, params_obs::Observable; dx=1)
     V[:, 1] .= V[:, end-1]
     V[:, end] .= V[:, 2]
 
-    # Update heatmap observable
-    return U
+    return U  # for visualization
 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/utils/constants.jl b/src/utils/constants.jl
index c06f46c..a07f4f0 100644
--- a/src/utils/constants.jl
+++ b/src/utils/constants.jl
@@ -1,6 +1,8 @@
 module Constants
 
-struct FHNParams
+abstract type PDEParams end
+
+struct FHNParams <: PDEParams
     N::Int
     dx::Float64 # grid spacing
     Du::Float64
@@ -15,7 +17,7 @@ struct FHNParams
         new(N, dx, Du, Dv, ϵ, a, b)
 end
 
-struct GSParams
+struct GSParams <: PDEParams
     N::Int       # grid size
     dx::Float64  # grid spacing
     Du::Float64  # diffusion rate U
@@ -25,6 +27,6 @@ struct GSParams
 
 end
 
-export FHNParams, GSParams
+export PDEParams, FHNParams, GSParams
 
 end # module Constants
diff --git a/src/utils/laplacian.jl b/src/utils/laplacian.jl
index e04d5a2..814573a 100644
--- a/src/utils/laplacian.jl
+++ b/src/utils/laplacian.jl
@@ -1,36 +1,16 @@
 module Laplacian
 
-"""
-    laplacian(U::Matrix{Float64}, N::Int, h::Float64)
-
-    Computes the discrete 2D Laplacian of a matrix `U` using a 5-point stencil
-    and circular boundary conditions.
-
-    # Arguments
-    - `U::Matrix{Float64}`: The input 2D matrix representing the field or image.
-    - `N::Int`: Integer
-    - `h::Float64`: The spatial step size or grid spacing between points in the discretization.
-
-    # Returns
-    - `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
-
-function laplacian5(f, dx)
+function laplacian(U::AbstractMatrix{<:Real}, dx::Real)
     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
+    center = U[2:end-1, 2:end-1]
+    up     = U[1:end-2, 2:end-1]
+    down   = U[3:end,   2:end-1]
+    left   = U[2:end-1, 1:end-2]
+    right  = U[2:end-1, 3:end]
+
+    return (up .+ down .+ left .+ right .- 4 .* center) ./ h2
 end
 
 export laplacian, laplacian5
 
-end # module Laplacian
\ No newline at end of file
+end # module Laplacian
diff --git a/src/utils/templates.jl b/src/utils/templates.jl
new file mode 100644
index 0000000..7298bf4
--- /dev/null
+++ b/src/utils/templates.jl
@@ -0,0 +1,198 @@
+
+# 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=5e-6, Dv=2e-3, ϵ=0.025, a=0.6, b=0.15)
+    params = FHNParams(N=N, dx=1.0, Du=0.182, Dv=0.5, ϵ=0.05, a=0.2, b=2.0)
+
+    u0 = 0.05 .* (2 .* rand(N, N) .- 1)  # noise in [-0.05, 0.05]
+    v0 = 0.05 .* (2 .* rand(N, N) .- 1)
+
+    return params, vcat(vec(u0), vec(v0))
+end
+
+function coral_conditions(N)
+    # Turing-spot parameters
+    params = FHNParams(N=N, dx=1.0, Du=0.001, Dv=0.06, ϵ=0.05, a=0.0, b=1.2)
+
+    u0, v0 = coral_ic(N)
+
+    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
+
+function coral_ic(N)
+    u = fill(0.534522, N, N)
+    v = fill(0.381802, N, N)
+
+    for _ in 1:40  # place 15 noisy seeds
+        i, j = rand(10:N-10), rand(10:N-10)
+        u[i-2:i+2, j-2:j+2] .= -0.534522 .+ 0.2 * rand(5, 5)
+    end
+
+    return vec(u), vec(v)
+end
diff --git a/src/visualization.jl b/src/visualization.jl
index d2f8343..c5a6fb0 100644
--- a/src/visualization.jl
+++ b/src/visualization.jl
@@ -1,43 +1,8 @@
 module Visualization
 include("gray_scott_solver.jl")
+#include("fhn_solver.jl")
+
 using GLMakie, Observables, Makie
-using .GrayScottSolver
-"""
-    step_through_solution(sol::SolutionType, N::Int)
-
-    Function for visualization for the output of run_simulation
-
-    # Arguments
-    - `sol`: computed differential equation by run_simulation
-    - `N`: size of the N×N grid
-
-    # Returns
-    - ``: Displays created figure
-"""
-function step_through_solution(sol, N::Int)
-    fig = Figure(resolution=(600, 600))
-    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)
@@ -45,6 +10,20 @@ function coord_to_index(x, y, N)
     return ix, iy
 end
 
+"""
+    reset(U, V, heat_obs)
+
+    Resets heatmap to original state by replacing each cell.
+    Currently only places a square in the center
+
+    # Arguments
+    - `U`: Matrix filled with ones
+    - `V`: Empty matrix filled with zeros
+    - `heat_obs`: Heatmap observable
+
+    # Returns
+    - ``: resetted heatmap observable
+"""
 function reset!(U, V, heat_obs)
     U .= 1.0
     V .= 0.0
@@ -68,7 +47,7 @@ function param_box!(grid, row, labeltxt, observable::Observable)
     end
 end
 
-function build_ui(U, V, Du, Dv, F, k, params_obs, heat_obs)
+function build_ui(U, V, param_obs_map::NamedTuple, params_obs, heat_obs)
     reset!(U, V, heat_obs)
     fig = Figure(size=(800, 800))
     gh = GridLayout(fig[1, 1])
@@ -81,7 +60,7 @@ function build_ui(U, V, Du, Dv, F, k, params_obs, heat_obs)
     stepsize = Observable(30)
     spoint = select_point(ax.scene)
 
-    # # Controls
+    # 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)
@@ -95,10 +74,17 @@ function build_ui(U, V, Du, Dv, F, k, params_obs, heat_obs)
     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)
+    param_box!(textboxgrid, 1, "Du", param_obs_map.Du)
+    param_box!(textboxgrid, 2, "Dv", param_obs_map.Dv)
+
+    if haskey(param_obs_map, :F)
+        param_box!(textboxgrid, 3, "Feed", param_obs_map.F)
+        param_box!(textboxgrid, 4, "kill", param_obs_map.k)
+    elseif haskey(param_obs_map, :ϵ)
+        param_box!(textboxgrid, 3, "ϵ", param_obs_map.ϵ)
+        param_box!(textboxgrid, 4, "a", param_obs_map.a)
+        param_box!(textboxgrid, 5, "b", param_obs_map.b)
+    end
 
     # Timer and state for animation
     running = Observable(false)
@@ -159,6 +145,5 @@ function build_ui(U, V, Du, Dv, F, k, params_obs, heat_obs)
     return fig
 end
 
-
-export step_through_solution, build_ui
+export build_ui
 end