Differential Equations and Simulation

Simulating how a system changes over time — populations, chemicals, circuits, orbits — usually means solving ordinary differential equations (ODEs): du/dt = f(u, p, t).

A solver from scratch

The classic fourth-order Runge–Kutta method takes a few lines of generic Julia:

function rk4(f, u0, p, tspan, dt)
    t0, t1 = tspan
    ts = collect(t0:dt:t1)
    us = [u0]
    u = u0
    for t in ts[1:end-1]
        k1 = f(u, p, t)
        k2 = f(u .+ dt / 2 .* k1, p, t + dt / 2)
        k3 = f(u .+ dt / 2 .* k2, p, t + dt / 2)
        k4 = f(u .+ dt .* k3, p, t + dt)
        u = u .+ dt / 6 .* (k1 .+ 2 .* k2 .+ 2 .* k3 .+ k4)
        push!(us, u)
    end
    return ts, us
end

# exponential decay du/dt = -k u, exact solution exp(-k t)
decay(u, k, t) = -k * u
ts, us = rk4(decay, 1.0, 0.5, (0.0, 4.0), 0.1)
println(round(us[end], digits = 6), " vs exact ", round(exp(-0.5 * 4), digits = 6))

# Lotka–Volterra predator–prey model
function lotka(u, p, t)
    prey, pred = u
    α, β, δ, γ = p
    return [α * prey - β * prey * pred, δ * prey * pred - γ * pred]
end
ts, us = rk4(lotka, [10.0, 5.0], (1.1, 0.4, 0.1, 0.4), (0.0, 20.0), 0.01)
prey = first.(us)
println("prey range: ", round.(extrema(prey), digits = 2))

0.135335 vs exact 0.135335
prey range: (0.51, 13.68)

Because rk4 is generic, the same code works for a scalar, a vector, or any type that supports + and *.

DifferentialEquations.jl

For real work, use DifferentialEquations.jl (part of the SciML ecosystem): hundreds of solvers, automatic step-size control, stiff problems, events, stochastic and delay equations.

# illustrative — needs the DifferentialEquations package
using DifferentialEquations

function lotka!(du, u, p, t)
    prey, pred = u
    α, β, δ, γ = p
    du[1] = α * prey - β * prey * pred
    du[2] = δ * prey * pred - γ * pred
end

prob = ODEProblem(lotka!, [10.0, 5.0], (0.0, 20.0), (1.1, 0.4, 0.1, 0.4))
sol = solve(prob, Tsit5(), reltol = 1e-8)

sol(7.5)            # interpolated state at any time
sol.t[end], sol.u[end]

using Plots
plot(sol, labels = ["prey" "predator"])

• The in-place form f!(du, u, p, t) writes into du and avoids allocations.
• Tsit5() is a good default for non-stiff problems; Rodas5P() or QNDF() for stiff ones — or call solve(prob) and let it choose.
• sol is callable: sol(t) interpolates between steps.

The wider SciML ecosystem

• ModelingToolkit.jl — write models symbolically (equations and components), and it generates fast simulation code.
• SciMLSensitivity.jl — gradients of solutions with respect to parameters, for fitting and machine learning.
• Catalyst.jl — chemical reaction networks.
• DiffEqFlux.jl and neural ODEs — combine differential equations with neural networks ("scientific machine learning").