function ex19()
m = 75;
g = 9.81;
b = 0.1;
y0 = [5000; 0];

h = 10;
tE = 70;
f = @(t,y) [y(2); -g - (b/m)*y(2)^2*sign(y(2))];

t = (0:h:tE)';
N = length(t);
y = zeros(N,2);
y(1,:) = y0';

for I = 2:N
  y(I,:) = odeSchritt(f, t(I-1), y(I-1,:)', h);
end

[tM, yM] = ode45(f, t, y0);

% Ergebnismatrix zum schöneren Ausgeben 
OUT = [t, y(:,1), yM(:,1)];
fprintf('  t    x(Hand)   x(Matlab)\n');
fprintf('==========================\n');
fprintf('%3.0f   %8.2f   %8.2f\n', OUT');

end

%---------------------------------------------------------------
function yNeu = odeSchritt(f, t, y, h)
 k1 = f(t, y);
 k2 = f(t + h/3, y + (1/3)*h*k1);
 k3 = f(t + (2/3)*h, y + h*(-k1/3 + k2));
 k4 = f(t + h, y + h*(k1 - k2 + k3));
 yNeu = y + (h/8)*(k1 + 3*k2 + 3*k3 + k4);
end
