%
%	Execute a Runge Kutta method for the solution
% 	of an IVP with shape
%		y' = f(t, y); 	y(t0) = y0
%
% 	Arguments:
% 		A         (matrix)    Butcher tableau matrix.
% 		b         (vector)    Butcher tableau b.
% 		c         (vector)    Butcher tableau c.
% 		t0        (scalar)    Initial time.
% 		y0        (vector)    Initial value.
% 		f         (function)  Right-hand side of ODE.
% 		t_limit   (scalar)    Time limit to integrate up to.
% 		h         (scalar)    Step length (if equidistant).
%
% 	Returns:
% 		times       (vector)    Computed time steps.
% 		steps       (vector)    Computed solution steps.
%
function [times, steps] = run_rk(A, b, c, t0, y0, f, t_limit, h)
	% init output variables
	times = [t0];
	steps = [y0];

	% init iteration variables
	t = t0;
	y = y0;

	% perform steps
	while t < t_limit
		y = rk(A, b, c, t, y, h, f);
		if t + h >= t_limit
			h = t_limit - t;
			t = t_limit;
		else
			t = t+h;
		end
		% y = rk(A, b, c, t, y, h, f);
		times = [times, t];
		steps = [steps, y];
	end
end