Tutorial · Intermediate · 12 min

Solving differential equations

ode45 solves dy/dt = f(t, y) with an adaptive Runge–Kutta method. Learn single equations, systems and phase portraits.

1. One equation

   Write the right-hand side as a function of t and y. ode45 returns the time points and the solution.

   f = @(tt, y) -0.5 * y;
   [tt, y] = ode45(f, [0 10], 2);
   plot(tt, y, 'o-')
   title("y' = -0.5 y,  y(0) = 2")
   xlabel('t'), ylabel('y')

   Try it ↗

   Tip: The time variable t is reserved for animations, so these examples call it tt.

2. Check against the exact answer

   The exact solution is 2·e^−t/2. Plot both, and the error on a log scale.

   [tt, y] = ode45(@(tt, y) -0.5*y, [0 10], 2);
   exact = 2 * exp(-tt/2);
   subplot(2, 1, 1)
   plot(tt, y, 'o', tt, exact, '-')
   legend('ode45', 'exact')
   subplot(2, 1, 2)
   semilogy(tt, abs(y - exact) + eps)
   title('absolute error')

   Try it ↗

3. A system: the pendulum

   A second-order equation θ'' = −sin θ becomes a system for y = [θ; ω]. The function returns a column vector.

   f = @(tt, y) [y(2); -sin(y(1))];
   [tt, y] = ode45(f, [0 20], [2.5; 0]);
   plot(tt, y(:, 1), tt, y(:, 2), '--')
   legend('angle θ', 'angular velocity ω')
   xlabel('time')

   Try it ↗

4. Phase portrait

   Plot ω against θ for several starting angles, on top of the vector field.

   [TH, W] = meshgrid(linspace(-pi, pi, 21), linspace(-3, 3, 17));
   quiver(TH, W, W, -sin(TH))
   hold on
   for a = [0.5 1.5 2.5 3.1]
     [tt, y] = ode45(@(tt, y) [y(2); -sin(y(1))], [0 15], [a; 0]);
     plot(y(:, 1), y(:, 2), 'LineWidth', 1.5)
   end
   hold off
   xlabel('θ'), ylabel('ω'), title('Pendulum phase portrait')

   Try it ↗