%
%	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