Kuramoto model (ten lines)
Kuramoto model is a mathematical model of synchronization. Consider oscillators with random periods. (In tems of the following code, the random periods are denoted by 2\pi/\omega[i] and \omega[i]'s are generated by Normal(1, 2).) The model describes the synchronization of the oscillators:

If the strength K of synchronization is less than a certain critical value K_c, then no synchronization will occur.

Synchronization begins when the strength K of synchronization exceeds the ctitical value K_c.
The critical value K_c is determined by the probability distribution dist of \omega[i]'s: K_c = π/2/pdf(dist, mean(dist)).
The video below looks like fireflies are synchronizing their rhythmic flashes.
a = 1.3; tmax = 25.0; using Distributions, DifferentialEquations, Plots
kuramoto!(dθ, θ, param, t) = (N = length(θ); (K, ω) = param;
for i in 1:N dθ[i] = ω[i] + K*mean(sin(θ[j]  θ[i]) for j in 1:N) end)
(m, n) = (32, 16); (d, v) = (Normal(0, 2), 1.0); K_c = 2/π/pdf(d, 0)
θ₀ = 2π*rand(m, n); tspan = (0.0, tmax); K = a*K_c; ω = rand(d, m, n) .+ v
sol = solve(ODEProblem(kuramoto!, θ₀, tspan, (K, ω)))
anim = @animate for t in [fill(0, 20); 0:0.1:tmax; fill(tmax, 20)]
plot(size=(400, 220), title="Kuramoto model K = $(a)K_c: t = $t", titlefontsize=10)
heatmap!(sin.(sol(t))'; c=:bwr, colorbar=false, frame=false, ticks=false)
end; gif(anim, "kuramoto$(a).gif")
Result
For the readability, expand the above compressed code, in an appropriate environment, by
quote
a = 1.3; tmax = 25.0; using Distributions, DifferentialEquations, Plots
kuramoto!(dθ, θ, param, t) = (N = length(θ); (K, ω) = param;
for i in 1:N dθ[i] = ω[i] + K*mean(sin(θ[j]  θ[i]) for j in 1:N) end)
(m, n) = (32, 16); (d, v) = (Normal(0, 2), 1.0); K_c = 2/π/pdf(d, 0)
θ₀ = 2π*rand(m, n); tspan = (0.0, tmax); K = a*K_c; ω = rand(d, m, n) .+ v
sol = solve(ODEProblem(kuramoto!, θ₀, tspan, (K, ω)))
anim = @animate for t in [fill(0, 20); 0:0.1:tmax; fill(tmax, 20)]
plot(size=(400, 220), title="Kuramoto model K = $(a)K_c: t = $t", titlefontsize=10)
heatmap!(sin.(sol(t))'; c=:bwr, colorbar=false, frame=false, ticks=false)
end; gif(anim, "kuramoto$(a).gif")
end > Base.remove_linenums! > x > display("text/markdown", "```julia\n$x\n```")
Expanded Code
begin
a = 1.3
tmax = 25.0
using Distributions, DifferentialEquations, Plots
kuramoto!(dθ, θ, param, t) = begin
N = length(θ)
(K, ω) = param
for i = 1:N
dθ[i] = ω[i] + K * mean((sin(θ[j]  θ[i]) for j = 1:N))
end
end
(m, n) = (32, 16)
(d, v) = (Normal(0, 2), 1.0)
K_c = (2 / π) / pdf(d, 0)
θ₀ = (2π) * rand(m, n)
tspan = (0.0, tmax)
K = a * K_c
ω = rand(d, m, n) .+ v
sol = solve(ODEProblem(kuramoto!, θ₀, tspan, (K, ω)))
anim = @animate(for t = [fill(0, 20); 0:0.1:tmax; fill(tmax, 20)]
plot(size = (400, 220), title = "Kuramoto model K = $(a)K_c: t = $(t)", titlefontsize = 10)
heatmap!((sin.(sol(t)))'; c = :bwr, colorbar = false, frame = false, ticks = false)
end)
gif(anim, "kuramoto$(a).gif")
end