Write a MATLAB function [T, Y] = rk4(f, tspan, y0, n) that solves dy/dt = f(t, y) with initial condition y(t = a) = y0 on the interval [a, b] using n steps of the classical 4th order Runge-Kutta method.
Inputs:
f - a function f(t, y) that returns a column vector of the same length as y
tspan - a vector [a, b] with the initial and end values of t
y0 - a column vector of the initial values
n - number of steps to use
Outputs:
T - an (n+1) column vector containing the times
Y - an (n+1) by d matrix, where d is the length of y; Y(t, i) gives the i-th component of y at timestep t
For example, [T, Y] = rk4(@(t, y) -y, [0 1], 1, 2) should return:
T = [0
0.5
1]
Y = [1
0.60677
0.36817]
As another example, [T, Y] = rk4(@(t, y) [y(2); -y(1)], [0 1], [1; 0], 2) should return:
T = [0
0.5
1]
Y = [1 0
0.87760 -0.47917
0.54059 -0.84104]