%
% Execute a method for the solution of an IVP with shape
% y' = f(t, y); y(t0) = y0
%
% Arguments:
% t0 (scalar) Initial time.
% y0 (vector) Initial value.
% f (function) Right-hand side of ODE.
% nsteps (integer) Number of integration steps.
% h (scalar) Step length (if equidistant).
%
% Returns:
% timepoints (vector) Computed time steps.
% yapprox (vector) Computed solution approximations.
%
function [timepoints, yapprox] = blackbox_method(f, t0, y0, nsteps, h)
% init output variables
timepoints = zeros(1, nsteps+1);
timepoints(1) = t0;
yapprox = zeros(length(y0), nsteps+1);
yapprox(:,1) = y0;
% init iteration variables
t = t0;
y = y0;
% perform steps
for k = 1:nsteps
% intermediate computations for current step
k1 = f(t, y);
k2 = f(t + (1/5)*h, y + h*((1/5)*k1));
k3 = f(t + (3/10)*h, y + h*((3/40)*k1 + (9/40)*k2));
k4 = f(t + (4/5)*h, y + h*((44/45)*k1 - (56/15)*k2 + (32/9)*k3));
k5 = f(t + (8/9)*h, y + h*((19372/6561)*k1 - (25360/2187)*k2 + (64448/6561)*k3 - (212/729)*k4));
k6 = f(t + h, y + h*((9017/3168)*k1 - (355/33)*k2 + (46732/5247)*k3 + (49/176)*k4 - (5103/18656)*k5));
% perform step
y = y + h*((35/384)*k1 + (500/1113)*k3 + (125/192)*k4 - (2187/6784)*k5 + (11/84)*k6);
t = t + h;
timepoints(k+1) = t;
yapprox(:,k+1) = y;
end
end