How to Solve Differential Equations in MATLAB (ode45, ode15s, bvp4c)

Differential equation assignments usually boil down to three scenarios: standard initial value problems, stiff systems that crash normal solvers, and boundary value problems where conditions are split between two ends of a domain.

This guide walks through how to pick the right MATLAB solver, structure your equations into first-order systems, and write clean code that yields accurate plots.

Quick Solver Decision Guide

Part 1: Standard Initial Value Problems (ode45)

ode45 is the default solver in MATLAB. It uses a medium-order Runge-Kutta algorithm (4th/5th order) and works well for most smooth, non-stiff dynamic systems.

Problem Statement

Solve a second-order damped harmonic oscillator:

y'' + 0.5 y' + 9 y = 0, with y(0) = 2, y'(0) = 0, for t in [0, 10]

Step 1: Reduce to First-Order System

MATLAB solvers require a system of first-order differential equations. Define state variables:

y1 = y, y2 = dy/dt

Derivatives:

dy1/dt = y2
dy2/dt = -0.5 * y2 - 9 * y1

Step 2: MATLAB Code Implementation

function solve_oscillator()
    % Time span and initial conditions [y(0), y'(0)]
    tspan = [0, 10];
    y0 = [2; 0];

    % Solve using ode45
    [t, y] = ode45(@oscillator_system, tspan, y0);

    % Plot results
    figure;
    plot(t, y(:,1), 'b-', 'LineWidth', 1.5); hold on;
    plot(t, y(:,2), 'r--', 'LineWidth', 1.5);
    xlabel('Time (s)');
    ylabel('State Response');
    title('Damped Harmonic Oscillator (ode45)');
    legend('Position y(t)', 'Velocity y''(t)');
    grid on;
end

function dydt = oscillator_system(t, y)
    dydt = zeros(2, 1);
    dydt(1) = y(2);
    dydt(2) = -0.5*y(2) - 9*y(1);
end

Part 2: Stiff Systems (ode15s)

A system is stiff when it contains fast and slow time scales simultaneously. If you run ode45 on a stiff problem, time steps shrink to tiny numbers, causing execution to freeze. ode15s uses numerical differentiation formulas to solve stiff systems efficiently.

Problem Statement (Robertson Reaction Model)

dy1/dt = -0.04 * y1 + 1e4 * y2 * y3
dy2/dt = 0.04 * y1 - 1e4 * y2 * y3 - 3e7 * y2^2
dy3/dt = 3e7 * y2^2

Initial state: y(0) = [1, 0, 0] over t in [0, 4e5]

Step 1: MATLAB Code Implementation

function solve_stiff_system()
    tspan = [0, 4e5];
    y0 = [1; 0; 0];

    % ode15s is required due to extreme rate constant differences
    options = odeset('RelTol', 1e-4, 'AbsTol', [1e-6, 1e-10, 1e-6]);
    [t, y] = ode15s(@robertson_system, tspan, y0, options);

    % Plot concentration curves on log scale
    figure;
    semilogx(t, y(:,1), 'r-', 'LineWidth', 1.5); hold on;
    semilogx(t, y(:,2)*1e4, 'g--', 'LineWidth', 1.5);
    semilogx(t, y(:,3), 'b-.', 'LineWidth', 1.5);
    xlabel('Time (s, log scale)');
    ylabel('Concentration');
    title('Robertson Stiff Reaction Kinetics (ode15s)');
    legend('Species y1', 'Species y2 (x 10^4)', 'Species y3');
    grid on;
end

function dydt = robertson_system(t, y)
    dydt = zeros(3, 1);
    dydt(1) = -0.04*y(1) + 1e4*y(2)*y(3);
    dydt(2) = 0.04*y(1) - 1e4*y(2)*y(3) - 3e7*y(2)^2;
    dydt(3) = 3e7*y(2)^2;
end

Part 3: Boundary Value Problems (bvp4c)

Unlike initial value problems where all states are known at t = 0, boundary value problems give conditions at two different spatial points. bvp4c solves boundary value problems using finite difference algorithms.

Problem Statement (Heat Transfer in a Cooling Fin)

d2T/dx2 - 0.15 * (T - 25) = 0, for x in [0, 2]

Boundary Conditions: T(0) = 100°C and dT/dx(2) = 0

Step 1: MATLAB Code Implementation

function solve_cooling_fin()
    % Mesh initialization and initial guess
    xmesh = linspace(0, 2, 10);
    solinit = bvpinit(xmesh, @initial_guess);

    % Solve BVP
    sol = bvp4c(@fin_ode, @fin_bc, solinit);

    % Evaluate solution on fine grid
    x_fine = linspace(0, 2, 100);
    y_fine = deval(sol, x_fine);

    % Plot Temperature Profile
    figure;
    plot(x_fine, y_fine(1, :), 'k-', 'LineWidth', 1.5);
    xlabel('Fin Length x (m)');
    ylabel('Temperature T (°C)');
    title('Temperature Distribution along Cooling Fin (bvp4c)');
    grid on;
end

function dydx = fin_ode(x, y)
    dydx = [y(2); 0.15 * (y(1) - 25)];
end

function res = fin_bc(ya, yb)
    res = [ya(1) - 100; yb(2) - 0];
end

function yguess = initial_guess(x)
    yguess = [100; 0];
end

Summary Troubleshooting Table