Solving ODEs Numerically in MATLAB

Solving a single first-order ODE

General procedure:

1. Define a function representing the right-hand side of the ODE
2. Define the timespan and initial value
3. Call ode45
4. Plot the results

Example: Solve the ODE dy(t)/dt = α (1 + cos(t)) y(t) - γ y(t)² with the initial condition y(0) = y₀ over the time span 0 to 4π for α = 20, γ = 3, and y₀ = 20

Define a function representing the right-hand side of the ODE:

function yprime = rhsODE(t, y)
alpha = 20; gamma = 3;
yprime = alpha*(1 + cos(t))*y - gamma*y^2;

Solve the ODE and plot the results:

tspan = [0 4*pi];
y0 = 20;

[t_soln, y_soln] = ode45(@rhsODE, tspan, y0);

plot(t_soln, y_soln)

Solving a system of first-order ODEs

Solve the system of ODEs:

with the intial conditions x₁(0) = 4 and x₂(0) = 1 over the time span 0 to 8

Define a function representing the right-hand side of the ODE:

function zprime = systemODE(t, z)
x1 = z(1);
x2 = z(2);

x1prime = x1 - 2*x2;
x2prime =3*x1 - 2*x2;

zprime = [x1prime; x2prime];

Solve the ODE and plot the results:

tspan = [0 8];
z0 = [4, 1];

[t_soln, z_soln] = ode45(@systemODE, tspan, z0);

x1_soln = z_soln(:, 1);
x2_soln = z_soln(:, 2);
plot(t_soln, x1_soln, 'k:', t_soln, x2_soln, 'b')

Solving a higher-order ODE

Convert the higher-order ODE to standard form, i.e., a system of first-order ODEs.

Example: The motion of a damped spring-mass system can be described using a second-order ODE:

Convert to standard form:

In MATLAB:

function wprime = rhsODE(t, w)
c = ...
m = ...
k = ...

Example: The KdV equation is a third-order ODE:

Convert to standard form: