ODE solver and phase portraits

Write the system one equation per line: x' = y, y' = -x - 0.2*y. Higher-order equations are converted for you: θ'' = -sin(θ) - b*θ' creates the states θ and θ′. Any other name = value line is a parameter (parameters may use earlier ones, e.g. L = 2*g). Use t for time; Greek letters work as names. Functions: sin, cos, tan, exp, ln/log, sqrt, abs, sign, min, max, step, and more; ^ is a power.

Initial conditions look like θ = 2.5, θ' = 0. Pick a time span and a method:

Phase plane: choose the axes. The direction field and the nullclines (where each plotted component’s derivative is zero) are drawn with the other states held at their initial values. Click anywhere to launch a trajectory from that point, forward and optionally backward in time. Drag to pan, scroll to zoom, double-click to refit.

Equilibria (autonomous systems only) are found by Newton’s method from a grid of starting points across the visible range. The Jacobian is computed numerically and its eigenvalues (Faddeev–LeVerrier characteristic polynomial, Durand–Kerner roots) classify each point as a node, spiral, saddle, saddle-focus, or centre. A linear centre or a zero eigenvalue means the linearisation cannot decide stability.

3-D: for systems with three or more states (the Lorenz attractor, SIR), tick 3-D and drag to rotate.

Export the solution as CSV, or the phase portrait as PNG.

Type the rate dy/dt = g(t, y) using t and y (0.5*y, y*(1 - y), cos(t) - y), or pick one from Equations…. Each short segment shows the slope g(t, y) at that point: every solution curve runs along them.

Click the field to start from that point (t₀, y₀), or type t₀, y₀ and the end time T.

Δt (slider or box) is Euler’s step: yₙ₊₁ = yₙ + Δt·g(tₙ, yₙ). The orange path is Euler; the dashed curve is the true solution through the same point — the exact formula for the built-in equations (growth, decay, the stiff −3y, logistic, cooling, cos t), otherwise an adaptive RK45 reference curve, labelled as such. The readout gives Euler’s final value, the true one and the error; the table lists the first twelve steps as you would write them by hand.

Try dy/dt = −3y with Δt = 1: each step multiplies y by −2 and Euler explodes instead of decaying.

Euler — position and velocity both from the old state;

semi-implicit (symplectic) Euler — velocity first, then position with the new velocity;

velocity Verlet — the trapezoid-rule cousin used in molecular dynamics and games;

RK4 — classical fourth-order Runge–Kutta.

One-compartment PK (invented drug) — an oral dose absorbed from the gut at rate ka into one well-mixed compartment of volume V and eliminated at ke = ln 2 / t½. The drug and its numbers are invented; nothing here is a dose. The solution matches the Bateman equation C(t) = F·D·ka / (V(ka − ke)) · (e^(−ke·t) − e^(−ka·t)), peaking at t = ln(ka/ke)/(ka − ke). Repeated doses are not modelled (the system has no impulse input); change the initial amount to see proportionality.

FitzHugh–Nagumo action potential — a two-variable caricature of the action potential: a fast voltage with positive feedback and a slow recovery. Start a little above the threshold and you get a full spike; a little below and nothing — all-or-nothing. A steady current I = 0.5 makes it fire repeatedly: intensity is coded by frequency, not amplitude. Turn on the nullclines and click the phase plane to find the threshold. (FitzHugh 1961; Nagumo, Arimoto & Yoshizawa 1962.)

Minimal glucose model (toy) — shaped like Bergman’s minimal model (Bergman et al. 1979): glucose relaxes towards Gb on its own and faster as insulin action X builds; insulin rises while glucose is above h. The parameters are made up to show the shape; the units are there so the equations read sensibly, not to describe a person.