MATLAB Runge-Kutta 4th-order integration function

rk4.m

function y = rk4(f, y, t, dt, varargin)
%RK4 Runge-Kutta 4th order integration
%   RK4(F,Y,T,DT,...) integrates the differential equation y' = f(t,y) from
%   time T to time T+DT, where Y is the value of the solution vector at
%   time T. F is a function handle. For a scalar T and a vector Y, F(T,Y)
%   must return a column vector corresponding to f(t,y). Additional
%   arguments will be passed to the function F(T,Y,...).

k1 = f(t,        y,           varargin{:});
k2 = f(t + dt/2, y + k1*dt/2, varargin{:});
k3 = f(t + dt/2, y + k2*dt/2, varargin{:});
k4 = f(t + dt,   y + k3*dt,   varargin{:});

y = y + dt/6 * (k1 + 2*k2 + 2*k3 + k4);

end