function ex21()
t0 = 0;
t1 = 0.5;
tSpan = [t0 t1];
y0 = [1 0]';
epsrel = 0.1;

omega = 8.0;
f = @(t,y) [y(2); -omega^2*y(1)];

delta = 1e-3;

k1 = f(t0, y0);
r = norm(k1)/(max(norm(y0), delta));
h = 0.8/r * epsrel^(1/3);

t = t0;
y = y0;
tall = [t];
yall = [y'];     % y-Komponenten in einer Zeile

fprintf("   h      t      epslok\n");
while t < t1
  k2 = f(t + h/2, y + h/2*k1);
  k3 = f(t + 3/4*h, y + 3/4*h*k2);
  yneu = y + h/9*(2*k1 + 3*k2 + 4*k3);
  k4 = f(t + h, yneu);
  e = norm(h/72*(-5*k1 + 6*k2 + 8*k3 - 9*k4));
  
  yabs = min(norm(y), norm(yneu));
  epslok = e/max(yabs, delta);
  ok = "zu groß";
  if epslok <= epsrel
    ok = "ok";
    t = t + h;
    y = yneu;
    k1 = k4;
    tall = [tall; t];
    yall = [yall; y'];
  end
  
  fprintf("%.4f  %.4f  %f  %s\n", h, t, epslok, ok);
  h = h *min(5, 0.8 *(epsrel/epslok)^(1/3)); 
end

x = yall(:,1);
figure("Color", [1 1 1], "Position",[0 0 600 400]);
subplot(1,2,1);
plot(tall, x, "b-", tall, x, "k*");
xlim(tSpan);
title("Handcode","FontSize",12);

% Lösen  mit Matlab-Solver
options = odeset("RelTol", epsrel);
[tM, yM] = ode23(f, tSpan, y0, options);

x = yM(:,1);
subplot(1,2,2);
plot(tM, x, "b-", tM, x, "k*");
xlim(tSpan);
title("Matlab-Funktion","FontSize",12);

F = getframe(gcf);
imwrite(F.cdata, "bild66.png");

% Abweichung von der genauen Lösung
errHand = max(abs(cos(omega*tall) - yall(:,1)))
errMatlab = max(abs(cos(omega*tM) - yM(:,1)))

