function ex22()
mu = 10;
y0 = [1 0];
te = 60;
h = 0.01;
ts = 0:h:te;      % Vorgabe der Zwischenzeiten

f1 = @(t,y) f(t,y,mu);
J1 = @(t,y) J(t,y,mu);

% a. ohne Jacobimatrix
[ta, ya] = ode23t(f1, ts, y0);
xa = ya(:,1);

plot(ta, xa);
xlabel("t","FontSize",12);
ylabel("x","FontSize",12);
xlim([0 te])

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

% b. mit Jacobimatrix
options = odeset("Jacobian", J1);
[tb, yb] = ode23t(f1, ts, y0, options);
xb = yb(:,1);

% Abweichung a - b
diffAB = max(abs(xa - xb))

% c. mit eigenem Solver
[tc, yc] = odeTrapez(f1, [0 te], y0, h, J1);
xc = yc(:,1);

plot(ta, xa, tc, xc);
xlabel("t","FontSize",12);
ylabel("x","FontSize",12);
legend("ode23t", "odeTrapez", "Location", "southwest")
xlim([0 te])

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

% Abweichung a - c
diffAC = max(abs(xa - xc))

plot(ta, xc -xa);
xlabel("t","FontSize",12);
title("x_c - x_a","FontSize",12)
xlim([0 te])

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


%--------------------------------------------------------------------
function dy = f(t, y, mu)
% rechte Seite der van-der-Pool-Gleichung

dy1 = y(2);
dy2 = mu*(1 - y(1).^2).*y(2) - y(1);
dy = [dy1; dy2];

%--------------------------------------------------------------------
function df = J(t, y, mu)
% Jacobimatrix der van-der-Pool-Gleichung
df = [0, 1; -2*mu*y(1).*y(2) - 1, mu*(1 - y(1).^2)];
