Solving ODEs with DifferentialEquations.jl

using DifferentialEquations, Plots

# Logistic growth: du/dt = r*u*(1 - u/K)
function logistic!(du, u, p, t)
    r, K = p
    du[1] = r * u[1] * (1 - u[1] / K)
end

u0 = [1.0]                 # initial population
tspan = (0.0, 10.0)
p = (1.2, 100.0)           # r=1.2, K=100

prob = ODEProblem(logistic!, u0, tspan, p)
sol = solve(prob, Tsit5())  # adaptive 5th-order RK

# sol is callable & plotable
sol(0.0)            # [1.0]
sol(5.0)            # interpolated value at t=5
sol.u[end]          # final state

plot(sol, label="u(t)", xlabel="t", ylabel="population")

# Stiff system? Use a stiff solver
# sol_stiff = solve(prob, Rosenbrock23())

# Ensembles: simulate 100 perturbed runs in parallel
ensemble_prob = EnsembleProblem(prob, prob_func = (prob, i, repeat) ->
    remake(prob, u0 = [1.0 + 0.1rand()]))
sim = solve(ensemble_prob, Tsit5(), trajectories = 100)