diff --git a/animations/1.jl b/animations/1.jl
deleted file mode 100644
index e07ce29..0000000
--- a/animations/1.jl
+++ /dev/null
@@ -1,109 +0,0 @@
-using DrWatson
-@quickactivate "NonlinearDynamicsTextbook"
-include(srcdir("colorscheme.jl"))
-
-# %% The very first code snippet of the course
-using DynamicalSystems # load the library
-
-function lorenz_rule(u, p, t)
-    σ, ρ, β = p
-    x, y, z = u
-    dx = σ*(y - x)
-    dy = x*(ρ - z) - y
-    dz = x*y - β*z
-    return SVector(dx, dy, dz) # Static Vector
-end
-
-p  = [10.0, 28.0, 8/3] # parameters: σ, ρ, β
-u₀ = [0, 10.0, 0]      # initial state
-# create an instance of a `DynamicalSystem`
-lorenz = ContinuousDynamicalSystem(lorenz_rule, u₀, p)
-
-T  = 100.0 # total time
-dt = 0.01  # sampling time
-A  = trajectory(lorenz, T; dt)
-
-# %% Animate evolution of trajectories in the Standard map
-using InteractiveDynamics
-using DynamicalSystems
-using GLMakie
-
-# Standard map trajectories
-ds = Systems.standardmap(; k = 1.0)
-u0s = [[θ, p] for θ ∈ 0:2π for p ∈ 0:2π]
-lims = ((0, 2π), (0, 2π))
-
-figure, main, layout, obs = interactive_evolution(
-    ds, u0s; tail = 1000, lims,
-)
-main.xlabel = "θ"
-main.ylabel = "p"
-
-# %% Lorenz system trajectories
-using OrdinaryDiffEq: Tsit5, Vern9
-ds = Systems.lorenz()
-
-u0s =  [[10,20,40.0] .+ rand(3) for _ in 1:6]
-
-diffeq = (alg = Tsit5(), dtmax = 0.01)
-
-figure, main, layout, obs = interactive_evolution(
-    ds, u0s; tail = 1000, diffeq, colors = to_color.(COLORS)
-)
-
-# %% Add poincare plane
-import AbstractPlotting
-using AbstractPlotting.MakieLayout: LAxis
-o = AbstractPlotting.Point3f0(-25, 0, 0)
-w = AbstractPlotting.Point3f0(50, 0, 50)
-p = AbstractPlotting.FRect3D(o, w)
-
-# These lines are necessary for transparent planes
-a = AbstractPlotting.RGBAf0(0,0,0,0)
-c = AbstractPlotting.RGBAf0(0.2, 0.2, 0.8, 1.0)
-img = AbstractPlotting.ImagePattern([c a; a c]);
-AbstractPlotting.mesh!(main, p; color = img);
-
-# %% Plot Poincare sos
-psosplot = layout[:, 2] = LAxis(figure)
-psos = poincaresos(ds, (2, 0.0), 2000.0)
-AbstractPlotting.scatter!(psosplot, psos[:, 1], psos[:, 3])
-
-display(figure)
-
-# %% Henon heiles system trajectories
-ds = Systems.henonheiles()
-
-u0s = [[0.0, -0.25, 0.42081, 0.0],
-[0.0, 0.1, 0.5, 0.0],
-[0.0, -0.31596, 0.354461, 0.0591255]]
-
-diffeq = (alg = Vern9(), dtmax = 0.01)
-idxs = (1, 2)
-
-figure, main, layout, obs = interactive_evolution(
-    ds, u0s; idxs, tail = 2000, diffeq, colors = COLORS,
-)
-main.figure[AbstractPlotting.Axis][:names, :axisnames] = ("q₁", "q₂", "p₂")
-
-# %% Poincare brainscanning application
-using GLMakie, DynamicalSystems, InteractiveDynamics
-using OrdinaryDiffEq
-
-ds = Systems.henonheiles()
-diffeq = (alg = Vern9(),)
-u0s = [
-    [0.0, -0.25, 0.42081, 0.0],
-    [0.0, 0.1, 0.5, 0.0],
-    [0.0, -0.31596, 0.354461, 0.0591255]
-]
-trs = [trajectory(ds, 10000, u0; diffeq...)[:, SVector(1,2,3)] for u0 ∈ u0s]
-for i in 2:length(u0s)
-    append!(trs[1], trs[i])
-end
-
-# Inputs:
-j = 1 # the dimension of the plane
-tr = trs[1]
-
-brainscan_poincaresos(tr, j)
diff --git a/animations/11.jl b/animations/11.jl
deleted file mode 100644
index 0b249ac..0000000
--- a/animations/11.jl
+++ /dev/null
@@ -1,53 +0,0 @@
-using DrWatson
-@quickactivate "NonlinearDynamicsTextbook"
-include(srcdir("style.jl"))
-
-using DynamicalSystems, PyPlot
-
-# The solution explicitly assumes δL = 1
-∂²(u, i, L) = u[mod1(i-1, L)] -2u[i] + u[mod1(i+1, L)]
-function ∂⁴(u, i, L)
-    u[mod1(i-2, L)] +
-    -4*u[mod1(i-1, L)] +
-    6u[i] +
-    -4u[mod1(i+1, L)] +
-    u[mod1(i+2, L)]
-end
-
-function swift_hohenberg(T, r, dt, u0)
-    L = length(u0)
-    U = zeros(T, length(u0))
-    U[1, :] .= u0
-    for j in 1:T-1
-        @show j
-        for i in 1:L
-            U[j+1, i] = U[j, i] + dt*(
-                        (r-1)*U[j, i] -
-                        2∂²(view(U, j, :), i, L) -
-                        ∂⁴(view(U, j, :), i, L) - U[j, i]^3
-            )
-        end
-    end
-    return U
-end
-
-dt = 0.05
-u0 = rand(100)
-r = 0.5
-T = 1000
-
-U = swift_hohenberg(T, r, dt, u0)
-
-figure()
-imshow(U')
-xlabel("time")
-ylabel("space")
-gca().set_aspect("auto")
-colorbar()
-tight_layout(;pad = 0.25)
-
-# figure()
-# plot(U[end, :])
-# xlabel("space")
-# ylabel("u")
-# tight_layout(;pad = 0.25)
diff --git a/animations/11/Brusselator.jl b/animations/11/Brusselator.jl
new file mode 100644
index 0000000..ef7cfed
--- /dev/null
+++ b/animations/11/Brusselator.jl
@@ -0,0 +1,38 @@
+# This file needs `src/brusselator_produce.jl` to have run with appropriate parameters.
+# The `saveat` parameter of that file is the time output of the frames.
+using DrWatson
+@quickactivate "NonlinearDynamicsTextbook"
+include(srcdir("colorscheme.jl"))
+using Random, OrdinaryDiffEq
+import GLMakie
+
+a = 9.0
+b = 10.2
+L = 50.0 # size of domain
+N = 250  # no. of grid points
+d = 1.7
+
+prefix = "brusselator2D_dense"
+fname = savename(prefix, @strdict(a, b, d, L, N, ic), "jld2")
+data = load(datadir("Brusselator", fname))
+uout = data["u"]
+tvec = data["t"]
+
+heatobs = GLMakie.Observable(uout[1][1,:,:])
+tobs = GLMakie.Observable(0.0)
+titobs = GLMakie.lift(t -> "t = $(t)", tobs)
+
+fig = GLMakie.Figure(resolution = (600, 550))
+ax = fig[1,1] = GLMakie.Axis(fig; title = titobs)
+hmap = GLMakie.heatmap!(ax, heatobs; colormap = :inferno, colorrange = (0, 1.5))
+cb = GLMakie.Colorbar(fig[1, 2], hmap; width = 20)
+display(fig)
+
+GLMakie.record(
+        fig, projectdir("animations", "11", "brusselator_d=$(d)_b=$(b).mp4"), 
+        1:length(tvec); framerate = 15
+    ) do i
+    
+    tobs[] = tvec[i]
+    heatobs[] = uout[i][1,:,:]
+end
diff --git a/animations/11/brusselator_d=1.7_b=10.2.mp4 b/animations/11/brusselator_d=1.7_b=10.2.mp4
new file mode 100644
index 0000000..e2349f4
Binary files /dev/null and b/animations/11/brusselator_d=1.7_b=10.2.mp4 differ
diff --git a/animations/11/brusselator_d=2.0_b=9.8.mp4 b/animations/11/brusselator_d=2.0_b=9.8.mp4
new file mode 100644
index 0000000..22ada9f
Binary files /dev/null and b/animations/11/brusselator_d=2.0_b=9.8.mp4 differ
diff --git a/animations/11/brusselator_d=3.0_b=10.2.mp4 b/animations/11/brusselator_d=3.0_b=10.2.mp4
new file mode 100644
index 0000000..2afdbd8
Binary files /dev/null and b/animations/11/brusselator_d=3.0_b=10.2.mp4 differ
diff --git a/animations/11/brusselator_d=3.0_b=9.8.mp4 b/animations/11/brusselator_d=3.0_b=9.8.mp4
new file mode 100644
index 0000000..588dc7c
Binary files /dev/null and b/animations/11/brusselator_d=3.0_b=9.8.mp4 differ
diff --git a/animations/11/brusselator_d=4.0.mp4 b/animations/11/brusselator_d=4.0.mp4
new file mode 100644
index 0000000..255d72e
Binary files /dev/null and b/animations/11/brusselator_d=4.0.mp4 differ
diff --git a/animations/11/brusselator_d=4.0_b=9.8.mp4 b/animations/11/brusselator_d=4.0_b=9.8.mp4
new file mode 100644
index 0000000..78aa7eb
Binary files /dev/null and b/animations/11/brusselator_d=4.0_b=9.8.mp4 differ
diff --git a/animations/11/fitzhigh_various.jl b/animations/11/fitzhigh_various.jl
new file mode 100644
index 0000000..80a9e88
--- /dev/null
+++ b/animations/11/fitzhigh_various.jl
@@ -0,0 +1,139 @@
+using DrWatson
+@quickactivate "NonlinearDynamicsTextbook"
+include(srcdir("colorscheme.jl"))
+using Random, OrdinaryDiffEq
+import GLMakie
+
+function FitzHugh_Nagumo_ODE(du,u,p,t)
+    a, b, d, epsilon, hsq6 = p
+    uu = @view u[1,:,:]
+    vv = @view u[2,:,:]
+
+    # This code can be optimized a lot if time allows
+    # Create padded u matrix to incorporate Newman boundary conditions - 9-point stencil
+    uuu=[ [uu[2,2] uu[2,:]' uu[2,end-1] ] ; [ uu[:,2] uu uu[:,end-1] ] ; [ uu[end-1,2] uu[end-1,:]' uu[end-1,end-1] ] ]
+    diff_term = d .* ( 
+        4 .* (uuu[2:end-1,1:end-2] .+ uuu[2:end-1,3:end] .+ 
+            uuu[3:end,2:end-1] .+ uuu[1:end-2,2:end-1] ) .+ 
+            uuu[3:end,3:end] .+ uuu[3:end,1:end-2] .+ uuu[1:end-2,3:end] .+ 
+            uuu[1:end-2,1:end-2] .- 20 .*uu  
+    ) ./ hsq6  
+
+    du[1,:,:] = a .* uu .*(1 .-uu).*(uu.-b) .- vv .+ diff_term
+    du[2,:,:] = epsilon .* (uu .- vv )
+end
+
+# --------------------------  main program  -----------------------------
+
+@time begin
+
+# parameters
+L = 300   # size of domain
+N = 150   # no. of grid points
+h = L/N # spatial steps size
+hsq6 = 6*h*h
+
+# FHN parameters
+a = 3.    
+b = 0.2   
+d = 1.   #  diffusion
+epsilon = 0.01 
+
+# open figure
+kplt = 0  # no. of subplot
+
+allsols = []
+allts = []
+
+for ic = 1:2    # initial conditions ( = rows of the plot )
+
+println("ic=",ic)
+
+uu = zeros(N,N)
+vv = zeros(N,N)   
+Random.seed!(13371)
+uu = uu .+ 0.01*randn(N,N)    # random perturbation
+
+# local stimulius = concentric waves
+if ic == 1 
+   uu[40:41,40:41] .= 0.999
+   uu[75:76,75:76] .= 0.999
+end
+
+# local stimulius = spiral waves
+if ic == 2
+   uu[1:end,10:11] .= 0.999
+end   
+
+# initial values   
+u0 = zeros(2,N,N)
+u0[1,:,:] = uu
+u0[2,:,:] = vv 
+
+p = a, b, d, epsilon, hsq6 
+
+if ic == 1
+    saveat = 0.0:2:600
+    tspan = (saveat[1], saveat[end])
+
+    prob = ODEProblem(FitzHugh_Nagumo_ODE, u0, tspan, p)
+    sol = solve(prob, Tsit5(); reltol=1e-6, abstol=1e-6, maxiters = 1e7, saveat)
+    tvec = sol.t
+    uout = [u[1, :, :] for u in sol.u]
+    push!(allsols, uout)
+    push!(allts, tvec)
+
+elseif ic == 2
+    saveat = 0.0:2.0:300.0
+    tspan = (saveat[1], saveat[end])
+    prob = ODEProblem(FitzHugh_Nagumo_ODE, u0, tspan, p)
+    sol = solve(prob, Tsit5(); reltol=1e-6, abstol=1e-6, maxiters = 1e7, saveat)
+
+    uout = [u[1, :, :] for u in sol.u]
+    # second perturbation line
+    u = sol.u[end]
+    u[1, 55:73,1:70] .= 0.999
+    saveat = 300.0:2:900
+    tspan = (saveat[1], saveat[end])
+    prob = ODEProblem(FitzHugh_Nagumo_ODE, u, tspan, p)
+    sol = solve(prob, Tsit5(); reltol=1e-6, abstol=1e-6, maxiters = 1e7, saveat)
+
+    tvec = sol.t
+    uout2 = [u[1, :, :] for u in sol.u]
+    append!(uout, uout2)
+    push!(allsols, uout)
+    push!(allts, tvec)
+
+end
+end # ic loop
+end # cputime
+
+
+# %% Animate
+names = ("concentric", "plane_to_spiral")
+
+for i in (1,2)
+    name = names[i]
+    uout = allsols[i]
+    tvec = allts[i]
+
+    heatobs = GLMakie.Observable(uout[1])
+    tobs = GLMakie.Observable(0.0)
+    titobs = GLMakie.lift(t -> "t = $(t)", tobs)
+
+    fig = GLMakie.Figure(resolution = (600, 550))
+    ax = fig[1,1] = GLMakie.Axis(fig; title = titobs)
+    hmap = GLMakie.heatmap!(ax, heatobs; colormap = :tokyo, colorrange = (-0.2, 1))
+    cb = GLMakie.Colorbar(fig[1, 2], hmap; width = 20)
+    display(fig)
+
+    GLMakie.record(
+        fig, projectdir("animations", "11", "fitzhugh_$(name).mp4"), 
+        1:length(tvec); framerate = 15
+    ) do i
+    
+        tobs[] = tvec[i]
+        heatobs[] = uout[i]
+    end
+
+end
diff --git a/animations/11/fitzhugh_concentric.mp4 b/animations/11/fitzhugh_concentric.mp4
new file mode 100644
index 0000000..9876d38
Binary files /dev/null and b/animations/11/fitzhugh_concentric.mp4 differ
diff --git a/animations/11/fitzhugh_plane_to_spiral.mp4 b/animations/11/fitzhugh_plane_to_spiral.mp4
new file mode 100644
index 0000000..f9fa38e
Binary files /dev/null and b/animations/11/fitzhugh_plane_to_spiral.mp4 differ
diff --git a/animations/11/fitzhugh_spiral.jl b/animations/11/fitzhugh_spiral.jl
new file mode 100644
index 0000000..50c13b6
--- /dev/null
+++ b/animations/11/fitzhugh_spiral.jl
@@ -0,0 +1,28 @@
+# uses file produced by `11.6_produce.jl`.
+using DrWatson
+@quickactivate "NonlinearDynamicsTextbook"
+include(srcdir("colorscheme.jl"))
+using Random, OrdinaryDiffEq
+import GLMakie
+
+data = wload(datadir("FitzHugh", "spiralwave.jld2"))
+@unpack uout, tvec, params = data
+
+heatobs = GLMakie.Observable(uout[1][1,:,:])
+tobs = GLMakie.Observable(0.0)
+titobs = GLMakie.lift(t -> "t = $(t)", tobs)
+
+fig = GLMakie.Figure(resolution = (600, 550))
+ax = fig[1,1] = GLMakie.Axis(fig; title = titobs)
+hmap = GLMakie.heatmap!(ax, heatobs; colormap = :tokyo, colorrange = (-0.2, 1))
+cb = GLMakie.Colorbar(fig[1, 2], hmap; width = 20)
+display(fig)
+
+GLMakie.record(
+        fig, projectdir("animations", "11", "fitzhugh_spiralwave.mp4"), 
+        1:length(tvec); framerate = 15
+    ) do i
+    
+    tobs[] = tvec[i]
+    heatobs[] = uout[i][1,:,:]
+end
diff --git a/animations/11/fitzhugh_spiralwave.mp4 b/animations/11/fitzhugh_spiralwave.mp4
new file mode 100644
index 0000000..1ccafaf
Binary files /dev/null and b/animations/11/fitzhugh_spiralwave.mp4 differ
diff --git a/animations/11/scroll_wave.mp4 b/animations/11/scroll_wave.mp4
new file mode 100644
index 0000000..347d291
Binary files /dev/null and b/animations/11/scroll_wave.mp4 differ
diff --git a/animations/2/quasiperiodic.jl b/animations/2/quasiperiodic.jl
index 2d952b4..841df39 100644
--- a/animations/2/quasiperiodic.jl
+++ b/animations/2/quasiperiodic.jl
@@ -2,8 +2,8 @@
 using DrWatson
 @quickactivate "NonlinearDynamicsTextbook"
 include(srcdir("colorscheme.jl"))
-using GLMakie, DynamicalSystems, InteractiveDynamics
-
+using DynamicalSystems, InteractiveDynamics
+import GLMakie
 using OrdinaryDiffEq
 
 R = 2.0
@@ -30,18 +30,18 @@ end
 ds = ContinuousDynamicalSystem(quasiperiodic_f, [0.0, 0.0, 0.0], nothing)
 frequencies² = [7, 9]
 u0s = [[0, 0, √x] for x ∈ frequencies²]
-diffeq = (alg = Tsit5(), dtmax = 0.01)
+diffeq = (alg = Tsit5(), adaptive = false, dt = 0.01)
 
 lims = ((-R-r, R+r), (-R-r, R+r), (-3r, 3r))
 
-plotkwargs = [(linewidth = i^2*1.0,) for i in 1:length(u0s)]
+plotkwargs = [(linewidth = i^2*2.0,) for i in 1:length(u0s)]
 
 fig, obs = interactive_evolution(
     ds, u0s; tail = 20000, diffeq, colors = COLORSCHEME[1:length(u0s)], 
     transform = torus, lims, plotkwargs
 )
 
-main = content(fig[1,1])
+main = GLMakie.content(fig[1,1])
 
 freqs = collect("√"*string(s) for s in frequencies²)
 main.title = "Frequency ratios: "*join(freqs, ", ")
@@ -58,8 +58,14 @@ V = range(-π, π; length = 50)
 x1 = [(R + r*cos(θ))*cos(φ)      for θ in U, φ in V]
 y1 = [(R + r*cos(θ))*sin(φ)      for θ in U, φ in V]
 z1 = [r*sin(θ)                   for θ in U, φ in V]
-wireframe!(main, x1,y1,z1, shading = false, color = COLORS[3], linewidth = 0.1)
+GLMakie.wireframe!(main, x1,y1,z1;
+    shading = false, color = BLACK, linewidth = 0.1
+)
 
 
 # %%
 record_interaction(fig, string(@__FILE__)[1:end-2]*"mp4"; total_time = 10)
+
+# %% Clean figure for torus in book
+GLMakie.hidedecorations!(main)
+GLMakie.hidespines!(main)
\ No newline at end of file
diff --git a/animations/3.jl b/animations/3.jl
deleted file mode 100644
index 269aa06..0000000
--- a/animations/3.jl
+++ /dev/null
@@ -1,110 +0,0 @@
-using DrWatson
-@quickactivate "NonlinearDynamicsTextbook"
-include(srcdir("style.jl"))
-using DynamicalSystems, InteractiveDynamics
-import GLMakie
-using AbstractPlotting.MakieLayout
-
-# %% Sensitive dependence demonstration for Lorenz-63
-using OrdinaryDiffEq: Tsit5, Vern9
-ds = Systems.lorenz()
-u0 = 3ones(dimension(ds))
-u0s =  [u0 .+ 1e-3i for i in 1:3]
-
-diffeq = (alg = Vern9(), dtmax = 0.01)
-
-figure, main, layout, obs = interactive_evolution(
-    ds, u0s; tail = 100, diffeq, colors = COLORS,
-)
-
-
-# %% stretching and folding for logistic
-lo = Systems.logistic()
-X = trajectory(lo, 1000, Ttr = 100)
-# close("all")
-fig = figure()
-# ax = subplot2grid((1,3), (0,0))
-ax = subplot(121)
-ax.scatter(X[1:end-1], X[2:end], color = COLORS[1], s = 5, label = "\$x_{i+2}\$")
-ax.set_xlabel("\$x_{i}\$")
-ax.set_yticks([0, 0.5, 1])
-ax.scatter(X[1:end-2], X[3:end], color = COLORS[2], s = 5, label = "\$x_{i+2}\$")
-ax.legend(markerscale = 5)
-
-using3D()
-# ax2 = subplot2grid((1,3), (0,1), projection = "3d", colspan = 2)
-ax2 = subplot(122, projection = "3d")
-ax2.scatter(X[1:end-2], X[2:end-1], X[3:end], color = COLORS[1], s = 5)
-
-ax2.set_xticklabels([])
-ax2.set_yticklabels([])
-ax2.set_zticklabels([])
-ax2.set_xlabel("\$x_{i}\$", labelpad = -10)
-ax2.set_ylabel("\$x_{i+1}\$", labelpad = -10)
-ax2.set_zlabel("\$x_{i+2}\$", labelpad = -10)
-
-# fold arrow
-s = (0.15, 0.5, 1.2)
-e = (0.85, 0.5, 1.2)
-ax2.quiver3D(e..., (s .- e)..., color = COLORS[3])
-ax2.text3D(0.5, 0.1, 1.2, "Fold", color = COLORS[3])
-
-# stretch arrow
-s = (0.75, 0.35, 0.9)
-e = (0.95, 0.15, 0.1)
-ax2.quiver3D(e..., (s .- e)..., color = COLORS[5])
-
-s = (0.55, 0.45, 0.1)
-e = (0.75, 0.35, 0.9)
-ax2.quiver3D(e..., (s .- e)..., color = COLORS[5])
-ax2.text3D(1.1, 0.45, 0.1, "Stretch", color = COLORS[5])
-
-
-fig.tight_layout()
-fig.subplots_adjust(top = 0.95, wspace = 0.1)
-
-
-# %% 01-test for chaos for standard map
-using Statistics
-sm = Systems.standardmap(;k = 1.0)
-u1 = 0.1rand(2)
-u2 = u1 .+ [π, 0]
-
-fig, axs = subplots(2, 3;figsize = (figx, 2figy))
-c = 0.2
-for (i, u) in enumerate((u1, u2))
-    x = trajectory(sm, 10000, u)[:, 2]
-    axs[i, 1].plot(x[1:100]; color = "C$(i-1)", marker = "o", lw = 0.5)
-    p, q = ChaosTools.trigonometric_decomposition(x, c)
-    axs[i, 2].plot(p, q; color = "C$(i-1)", marker = "o", lw = 0.25, ms = 2)
-    Z = ChaosTools.mmsd(mean(x), p, q, length(x)÷10, c)
-    axs[i, 3].plot(Z; color = "C$(i-1)")
-end
-axs[1, 1].set_title("timeseries \$x_n\$")
-axs[1, 2].set_title("process \$(p_n, q_n)\$")
-axs[1, 3].set_title("sq. displ. \$Z_n\$")
-fig.tight_layout(;pad = 0.25)
-
-# %% 01-test chaos for henon heiles
-using Statistics
-sm = Systems.henonheiles()
-u0s = (
-    [0.0, -0.25, 0.42, 0.0], # chaotic
-    [0.0, 0.30266571044921875, 0.4205654433900762, 0.0], # periodic
-)
-
-fig, axs = subplots(2, 3;figsize = (figx, 2figy))
-c = 0.2
-dt = 1.0
-for (i, u) in enumerate(u0s)
-    x = trajectory(sm, 10000, u; dt)[:, 2]
-    axs[i, 1].plot(x[1:100]; color = "C$(i-1)", marker = "o", lw = 0.5)
-    p, q = ChaosTools.trigonometric_decomposition(x, c)
-    axs[i, 2].plot(p, q; color = "C$(i-1)", marker = "o", lw = 0.25, ms = 2)
-    Z = ChaosTools.mmsd(mean(x), p, q, length(x)÷10, c)
-    axs[i, 3].plot(Z; color = "C$(i-1)")
-end
-axs[1, 1].set_title("timeseries \$x_n\$")
-axs[1, 2].set_title("process \$(p_n, q_n)\$")
-axs[1, 3].set_title("sq. displ. \$Z_n\$")
-fig.tight_layout(;pad = 0.25)
diff --git a/animations/3/henon_folding.jl b/animations/3/henon_folding.jl
new file mode 100644
index 0000000..95214a2
--- /dev/null
+++ b/animations/3/henon_folding.jl
@@ -0,0 +1,60 @@
+using DrWatson
+@quickactivate "NonlinearDynamicsTextbook"
+include(srcdir("colorscheme.jl"))
+using InteractiveDynamics, DynamicalSystems
+import GLMakie
+
+xmax = 1.0
+xmin = -1.0
+ymin = -1.0
+ymax = 1.0
+d = 50
+const a = 1.4
+const b = 0.3
+
+ics = [
+    SVector(x, y) for x in range(xmin, xmax; length = d) for y in range(ymin, ymax; length = d)
+]
+
+f1(u) = SVector(u[1], 1 - a*u[1]^2 + u[2])
+f2(u) = SVector(b*u[1], u[2])
+f3(u) = SVector(u[2], u[1])
+function three_phases(ics)
+    p1 = f1.(ics)
+    p2 = f2.(p1)
+    p3 = f3.(p2)
+    return p1, p2, p3
+end
+
+function linear_transition!(o, p1, p2, steps, io)
+    for i in 1:steps # assumes we start with o[] = p1
+        o[] = p1 .+ (p2 .- p1) .* i/steps
+        GLMakie.recordframe!(io)
+        # sleep(0.01) # change to record frame
+    end
+end
+
+iterations = 10
+steps = 20
+
+fig = GLMakie.Figure(); display(fig)
+ax = GLMakie.Axis(fig[1,1])
+o = GLMakie.Observable(ics)
+GLMakie.scatter!(ax, o; 
+    color = COLORS[1], strokewidth = 0.5, strokecolor = :black, markersize = 5
+)
+ax.limits = ((-2.5,2.5),(-1.5,1.5))
+
+GLMakie.record(fig, string(@__FILE__)[1:end-2]*"mp4"; framerate = 30) do io
+    for j in 1:iterations
+        p1, p2, p3 = three_phases(ics)
+        linear_transition!(o, ics, p1, steps, io)
+        GLMakie.recordframe!(io)
+        linear_transition!(o, p1, p2, steps, io)
+        GLMakie.recordframe!(io)
+        linear_transition!(o, p2, p3, steps, io)
+        GLMakie.recordframe!(io)
+        global ics = p3
+        GLMakie.recordframe!(io)
+    end
+end
diff --git a/animations/3/henon_folding.mp4 b/animations/3/henon_folding.mp4
new file mode 100644
index 0000000..8b4b69c
Binary files /dev/null and b/animations/3/henon_folding.mp4 differ
diff --git a/animations/4.jl b/animations/4.jl
deleted file mode 100644
index eafc7f6..0000000
--- a/animations/4.jl
+++ /dev/null
@@ -1,85 +0,0 @@
-using DrWatson
-@quickactivate "NonlinearDynamicsTextbook"
-include(srcdir("colorscheme.jl"))
-using GLMakie, DynamicalSystems, InteractiveDynamics
-using Roots
-
-# %% Interactive bifurcations for 1D energy balance
-αtan(T) = 0.5 - (0.4/π)*atan((T-263)/2)
-dTdt(T, ε = 0.65, α=αtan, s= 1.0) = s*(1 - α(T)) - 1.6e-10 * ε * T^4
-
-Ts = 200:0.5:320.0
-
-function roots(ε)
-    f = T -> dTdt(T, ε)
-    r = Roots.find_zeros(f, Ts[1], Ts[end])
-    s = [to_color(COLORS[isodd(i) ? 4 : 3]) for i in 1:length(r)]
-    return [Point2f0(x, 0) for x in r], s # roots, stability
-end
-
-figure, layout = layoutscene(resolution = (1000, 800))
-
-ax = layout[1, :] = LAxis(figure)
-sll = labelslider!(figure, "ε =", 0:0.01:1.0; sliderkw = Dict(:startvalue => 0.65))
-layout[2, :] = sll.layout
-ε_observable = sll.slider.value
-
-# initialize plot
-fs = Observable(dTdt.(Ts))
-lines!(ax, Ts, fs; color = COLORS[2], linewidth = 4)
-lines!(ax, Ts, zero(Ts); color = COLORS[1])
-r, s = roots(0.65)
-rootvals = Observable(r)
-rootcols = Observable(s)
-scatter!(ax, rootvals; color = rootcols, markersize = 10)
-
-display(figure)
-
-on(ε_observable) do ε
-    fs[] = dTdt.(Ts, ε)
-    r, s = roots(ε)
-    rootvals[] = r
-    rootcols[] = s
-end
-
-# %% Interactive orbit diagram for logistic map
-using InteractiveDynamics, GLMakie
-using DynamicalSystems
-
-ds = Systems.logistic()
-p_min, p_max = 1.0, 4.0
-t = "orbit diagram for the logistic map"
-
-oddata = interactive_orbitdiagram(
-    ds, 1, p_min, p_max, 1;
-    parname = "r", title = t
-)
-
-# %% Cobweb for pomeaumannevile
-using InteractiveDynamics, GLMakie
-using DynamicalSystems
-
-using DynamicalSystems, PyPlot
-function pm_f(x, p, n)
-	γ, z = p
-	s =  -2(γ+1)
-    if x < -0.5
-		s*(x+0.5) - γ
-        # -4x - 3
-    elseif -0.5 ≤ x ≤ 0.5
-        @inbounds γ*x*(1 + abs(2x)^(z-1))
-    else
-        s*(x-0.5) + γ
-    end
-end
-
-rs = 0.7:0.01:1
-ds = DiscreteDynamicalSystem(pm_f, 0.3, [rs[1], 2.6])
-
-interactive_cobweb(
-    ds, rs, 1;
-    trajcolor = COLORS[1],
-    fkwargs = [(linewidth = 4.0, color = COLORS[i+1]) for i in 1:3],
-	xmin = -1,
-	xmax = 1
-)
diff --git a/animations/4/bifurcations_energybalance.jl b/animations/4/bifurcations_energybalance.jl
new file mode 100644
index 0000000..877ca16
--- /dev/null
+++ b/animations/4/bifurcations_energybalance.jl
@@ -0,0 +1,77 @@
+# %% Interactive bifurcations for 1D energy balance
+using DrWatson
+@quickactivate "NonlinearDynamicsTextbook"
+include(srcdir("colorscheme.jl"))
+using DynamicalSystems, InteractiveDynamics
+using Roots
+import GLMakie
+
+αtan(T) = 0.5 - (0.4/π)*atan((T-263)/2)
+dTdt(T, ε = 0.65, α=αtan, s= 1.0) = s*(1 - α(T)) - 1.6e-10 * ε * T^4
+
+Ts = 200:0.5:320.0
+εs = 0.2:0.01:1.0
+
+function roots(ε)
+    f = T -> dTdt(T, ε)
+    r = Roots.find_zeros(f, Ts[1], Ts[end])
+    s = [GLMakie.to_color(COLORS[isodd(i) ? 4 : 3]) for i in 1:length(r)]
+    return [Point2f0(x, 0) for x in r], s # roots, stability
+end
+
+fig = GLMakie.Figure(resolution = (800, 500))
+
+ax = fig[1, 1] = GLMakie.Axis(fig;
+    xlabel = GLMakie.L"T", ylabel = GLMakie.L"f(T)",
+    ylabelsize = 24, xlabelsize = 24
+)
+ax.limits = ((Ts[1], Ts[end]), (-0.3, 0.3))
+axb = fig[1, 2] = GLMakie.Axis(fig;
+    xlabel = GLMakie.L"\epsilon", ylabel = GLMakie.L"T^*",
+    ylabelsize = 24, xlabelsize = 24
+)
+sll = GLMakie.labelslider!(fig, "ϵ =", εs; sliderkw = Dict(:startvalue => 0.65))
+fig[2, :] = sll.layout
+ε_observable = sll.slider.value
+axb.limits = ((εs[1], εs[end]), (220, 300))
+
+
+# initialize plot
+fs = GLMakie.Observable(dTdt.(Ts))
+GLMakie.lines!(ax, Ts, fs; color = COLORS[2], linewidth = 4)
+GLMakie.hlines!(ax, 0; color = COLORS[1])
+r, s = roots(0.65)
+rootvals = GLMakie.Observable(r)
+rootcols = GLMakie.Observable(s)
+GLMakie.scatter!(ax, rootvals; color = rootcols, markersize = 10)
+
+stablefp = GLMakie.Observable(GLMakie.Point2f0[])
+unstablefp = GLMakie.Observable(GLMakie.Point2f0[])
+
+GLMakie.scatter!(axb, stablefp; color = COLORS[4])
+GLMakie.scatter!(axb, unstablefp; color = COLORS[3])
+
+display(fig)
+
+GLMakie.on(ε_observable) do ε
+    fs[] = dTdt.(Ts, ε)
+    r, s = roots(ε)
+    rootvals[] = r
+    rootcols[] = s
+
+    for (i, ρ) in enumerate(r)
+        p = GLMakie.Point2f0(ε, ρ[1])
+        if isodd(i)
+            push!(stablefp[], p)
+        else
+            push!(unstablefp[], p)
+        end
+    end
+    stablefp[] = stablefp[]
+    unstablefp[] = unstablefp[]
+end
+
+ε_observable[] = 0.65
+
+# %%
+record_interaction(fig, string(@__FILE__)[1:end-2]*"mp4"; total_time = 15, sleep_time = 2)
diff --git a/animations/4/bifurcations_energybalance.mp4 b/animations/4/bifurcations_energybalance.mp4
new file mode 100644
index 0000000..c919be5
Binary files /dev/null and b/animations/4/bifurcations_energybalance.mp4 differ
diff --git a/animations/5.jl b/animations/5.jl
deleted file mode 100644
index d756cb5..0000000
--- a/animations/5.jl
+++ /dev/null
@@ -1,2 +0,0 @@
-# TODO: I can adapt the code that puts boxes over henon and animate smaller and smaller
-# boxes, while adding elements to the right plot that are just scatter points
diff --git a/animations/9/coupled_roessler.jl b/animations/9/coupled_roessler.jl
new file mode 100644
index 0000000..0e46268
--- /dev/null
+++ b/animations/9/coupled_roessler.jl
@@ -0,0 +1,51 @@
+using DrWatson
+@quickactivate "NonlinearDynamicsTextbook"
+include(srcdir("colorscheme.jl"))
+
+using InteractiveDynamics
+using DynamicalSystems, OrdinaryDiffEq
+import GLMakie
+
+a = 0.2
+b = 0.2
+c = 5.7
+amu = 0.02
+ω1 = 1.1 + amu
+ω2 = 1.1 - amu
+k1 = 0.02
+k2 = k1
+
+y_init = [0.1, 0.2, 0., 0.11, 0.19, 0.1]
+
+diffeq = (alg = Tsit5(), adaptive = false, dt = 0.02, reltol = 1e-9, abstol = 1e-9)
+
+ds = Systems.coupled_roessler(y_init; ω1, ω2, a, b, c, k1, k2)
+
+ps = Dict(
+    1 => 1.0:0.01:1.5,
+    2 => 1.0:0.01:1.5,
+    5 => 0:0.01:6,
+    6 => 0:0.001:0.2,
+    7 => 0:0.001:0.2,
+)
+pnames = Dict(
+    1 => "ωx",
+    2 => "ωy",
+    5 => "c",
+    6 => "kx",
+    7 => "ky",
+)
+
+u0s =  [ds.u0]
+
+fig, obs = interactive_evolution_timeseries(
+    ds, u0s, ps; tail = 1000, diffeq, idxs = [1,4], pnames,
+    lims = ((-12,12), (-12,12))
+)
+
+ax = GLMakie.content(fig[1,1])
+ax.xlabel = "x1"
+ax.ylabel = "y1"
+
+# %%
+record_interaction(fig, string(@__FILE__)[1:end-2]*"mp4"; total_time = 30, sleep_time = 2)
\ No newline at end of file
diff --git a/animations/9/coupled_roessler.mp4 b/animations/9/coupled_roessler.mp4
new file mode 100644
index 0000000..e1c7e29
Binary files /dev/null and b/animations/9/coupled_roessler.mp4 differ
diff --git a/animations/9/generalized_sync.jl b/animations/9/generalized_sync.jl
new file mode 100644
index 0000000..bfaaf40
--- /dev/null
+++ b/animations/9/generalized_sync.jl
@@ -0,0 +1,61 @@
+using DrWatson
+@quickactivate "NonlinearDynamicsTextbook"
+include(srcdir("colorscheme.jl"))
+
+using InteractiveDynamics
+using DynamicalSystems, OrdinaryDiffEq
+import GLMakie
+
+function rossler_driving(u, p, t)
+    k = p[1]
+    x_1, x_2, x_3, y_1, y_2, y_3, z_1, z_2, z_3 = u
+
+    dx_1  =  - x_2 -x_3 
+    dx_2  =   x_1 + 0.2x_2  
+    dx_3  =  0.2 + x_3*(x_1-5.7)
+    (dx_1, dx_2, dx_3) = 6 .* (dx_1, dx_2, dx_3)
+    
+    dy_1  =  10*(- y_1 + y_2)
+    dy_2  =  28y_1 - y_2 - y_1*y_3 + k*x_2 
+    dy_3  =  y_1*y_2 - 2.666y_3
+        
+    dz_1  =  10*(- z_1 + z_2)
+    dz_2  =  28z_1 - z_2 - z_1*z_3 + k*x_2 
+    dz_3  =  z_1*z_2 - 2.666z_3
+
+    return SVector(dx_1, dx_2, dx_3, dy_1, dy_2, dy_3, dz_1, dz_2, dz_3)
+end
+
+x0 = [-1, -2, 0.1]
+y0 = [5.0, 5.0, 20.0]
+z0 = [-5.0, -5.0, 20.0]
+
+u0 = vcat(x0, y0, z0)
+k = 5.0
+
+ds = ContinuousDynamicalSystem(rossler_driving, u0, [k])
+
+ps = Dict(
+    1 => 0.0:0.01:12,
+)
+pnames = Dict(
+    1 => "k",
+)
+
+u0s =  [ds.u0]
+
+fig, obs = interactive_evolution_timeseries(
+    ds, u0s, ps; tail = 500, idxs = [4,7], pnames,
+    lims = ((-20,20), (-20,20))
+)
+
+GLMakie.Label(fig[:, :, GLMakie.Top()], "Roessler subsystem driving two Lorenz subsystems: y, z. Plotted are the first variables of each subsystem. See Eq.(9.9)", valign = :bottom, padding = (0, 0, 10, 0))
+GLMakie.content(fig[1,1]).xlabel = "y₁"
+GLMakie.content(fig[1,1]).ylabel = "z₁"
+GLMakie.content(fig[:, 2][1,1]).ylabel = "y₁"
+GLMakie.content(fig[:, 2][2,1]).ylabel = "z₁"
+GLMakie.content(fig[:, 2][2,1]).xlabel = "t"
+
+
+# %%
+record_interaction(fig, string(@__FILE__)[1:end-2]*"mp4"; total_time = 20, sleep_time = 2)
diff --git a/animations/9/generalized_sync.mp4 b/animations/9/generalized_sync.mp4
new file mode 100644
index 0000000..08dc5aa
Binary files /dev/null and b/animations/9/generalized_sync.mp4 differ
diff --git a/animations/9/kuramoto.jl b/animations/9/kuramoto.jl
new file mode 100644
index 0000000..b02895d
--- /dev/null
+++ b/animations/9/kuramoto.jl
@@ -0,0 +1,69 @@
+# %% Interactive bifurcations for 1D energy balance
+using DrWatson
+@quickactivate "NonlinearDynamicsTextbook"
+include(srcdir("colorscheme.jl"))
+using DynamicalSystems, OrdinaryDiffEq, Random
+using InteractiveDynamics: record_interaction
+import GLMakie
+
+D = 25
+omega_vec = randn(D)
+
+ds = Systems.kuramoto(D)
+diffeq = (alg = Tsit5(), adaptive = false, dt = 0.05)
+integ = integrator(ds; diffeq...)
+
+fig = GLMakie.Figure()
+display(fig)
+cmap = :curl
+axku = GLMakie.Axis(fig[1,1])
+
+
+# Other static elements
+GLMakie.hidedecorations!(axku)
+axku.aspect = GLMakie.DataAspect()
+GLMakie.lines!(axku, cos.(0:0.001:2π+0.002), sin.(0:0.001:2π+0.002);color = :black)
+
+# Plot balls
+phases = GLMakie.Observable(copy(integ.u))
+balls = GLMakie.lift(u -> [GLMakie.Point2f0(cos(φ), sin(φ)) for φ in u], phases)
+GLMakie.scatter!(axku, balls;
+    markersize  = 12, color = 1:D, colormap = cmap, strokewidth=2, strokecolor = :black
+)
+
+# Plot arrows
+arrows_end = GLMakie.lift(b -> 0.8b, balls)
+arrows_start = [GLMakie.Point2f0(0,0) for i in 1:D]
+GLMakie.arrows!(axku, arrows_start, arrows_end; 
+    color = 1:D, colormap = (cmap, 0.5), arrowsize = 10, linewidth = 4,
+)
+
+# Plot mean field Vector
+# add sliders and buttons
+fig[2, 1] = controllayout = GLMakie.GridLayout(tellwidth = false)
+run = controllayout[1, 1] = GLMakie.Button(fig; label = "run")
+Ks = 0.0:0.1:10.0
+kslider = GLMakie.labelslider!(fig, "K =", Ks; sliderkw = Dict(:startvalue => ds.p.K))
+controllayout[1, 2] = kslider.layout
+
+
+# run button functionality
+isrunning = GLMakie.Observable(false)
+GLMakie.on(run.clicks) do clicks; isrunning[] = !isrunning[]; end
+GLMakie.on(run.clicks) do clicks
+    @async while isrunning[]
+        GLMakie.isopen(fig.scene) || break # ensures computations stop if closed window
+        step!(integ)
+        phases[] = integ.u
+        yield()
+        sleep(0.001) # or `yield()` instead
+    end
+end
+
+GLMakie.on(kslider.slider.value) do val
+    integ.p.K = val
+    u_modified!(integ, true)
+end
+
+# %%
+record_interaction(fig, string(@__FILE__)[1:end-2]*"mp4"; total_time = 15, sleep_time = 2)
diff --git a/animations/9/kuramoto.mp4 b/animations/9/kuramoto.mp4
new file mode 100644
index 0000000..31c1090
Binary files /dev/null and b/animations/9/kuramoto.mp4 differ