function [ts, ys] = odeTrapez(f, tSpan, y0, h, J)
% Ode-Solver, der die implizite Trapezregel verwendet
% f         Function f(t, y) die die DGL definiert
% tSpan     Start- und Endzeit [t0 te]
% y0        Startwert y(t0) als Zeilenvektor (!)
% h         Zeitschritt (fest!)
% J         Jacobimatrix J(t, y) von f

N = length(y0);                        % Zahl der Variablen
ts = (tSpan(1):h:tSpan(2))';           % Speicher für die Ergebnisse besorgen
ys = zeros(length(ts), N);
ys(1,:) = y0;                          % Anfangswert ist erstes Ergebnis

idN = eye(N);
tol = 1e-5;                            % Genauigkeit des Newton-Solvers

for I = 1:(length(ts)-1)
  t = ts(I);
  y = ys(I,:)';
  ts(I+1) = t + h;
  
  F = @(x) x - (h/2)*f(t+h, x) - y - (h/2)*f(t,y); % Hilfsfunktion für Newton
  DF = @(x) idN - (h/2)*J(t+h, x);                 % Jacobimatrix dazu
  ys(I+1,:) = solveNewton(F, DF, y, tol);          % Startwert: aktuelles y
end
