MATLAB TUTORIAL for the Second Course, part 2.5: Runge--Kutta of order 4

Runge--Kutta of order 4

To approximate the solution of the system of differential equations

\[ \begin{split} \dot{x}_1 (t) &= f_1 \left( t, x_1 (t) , \ldots , x_n (t) \right) , \\ \dot{x}_2 (t) &= f_2 \left( t, x_1 (t) , \ldots , x_n (t) \right) , \\ \vdots & \qquad \vdots \\ \dot{x}_n (t) &= f_n \left( t, x_1 (t) , \ldots , x_n (t) \right) , \end{split} \]

with \( x_1 (a) =\alpha_1 , x_2 (a) = \alpha_2 , \ldots , x_n (a) = \alpha_n \) over the interval [a,b].

function [t,z]=rks4(F,a,b,za,m)
% Input:     F is the system input as a string 'F'
%            a and b are the end points of the interval
%            za=[x(a) y(a)] are the initial conditions
%            m is the number of steps
% Output:    t is the vector of steps
%            z=[x1(t) ... xn(t)] where xk(t) is the approximation
%             to the k-th dependent variable
h=(b-a)/m;
t=zeros(1,m+1);
z=zeros(m+1,length(za));
t=a:h:b;
z(1,:)=za;
for j=1:m
    k1=h*feval(F,t(j),z(j,:));
    k2=h*feval(F,t(j)+h/2,z(j,:)+k1/2);
    k3=h*feval(F,t(j)+h/2,z(j,:)+k2/2);
    k4=h*feval(F,t(j)+h,z(j,:)+k3);
    z(j+1,:)=z(j,:)+(k1+2*k2+2*k3+k4)/6;
end