function [y,t] = rk4(u0,T,H)
g = 20.0; %tyngdaccelerationen
alpha = pi/2;
F = g/(1 + H)^2; % drivkraft
y = []; % lista för uträknade värden av f
u = u0;
h = 0.01; % steglängd
y = [y, u];

f = @(u) [u(2),F*cos(alpha) - g/u(1)^2 + u(1)*((u(4)^2)), u(4), (F*sin(alpha) - 2*(u(2)*u(4)))/u(1)]';

%rk4
for t=h:h:T
k1 = f(u);
k2 = f(u + k1.*(h/2));
k3 = f(u + k2.*(h/2));
k4 = f(u + k3.*h);

u = u + (k1 + k2.*2 + k3.*2 + k4).*(h/6);
y = [y, u];
end
y = y';
t = 0:h:T;
t = t';