function [t,y,tsn] = imEuler(of, ti, ts, y0)

% impliziter Euler (Rueckwaerts-Euler)
% of Function Handle der Differentialgleichung
% ti Berechnungsintervall
% ts Schrittweite
% y0 Anfangsbedingung
% t,y Rueckgabewerte Funktion
% tsn = korrigierte Schrittweite

% Hinweis: Parameter an die Funktion of koennte mittels varargin uebergeben werden  

% W. Schweizer
% Simulation physikalischer Systeme
% 05/2016
% vers. 1.0

% Anpassen der korrekte Schrittweite
dt = ti(2)-ti(1);
n = round(dt/ts);   % Anzahl der Berechnungsschritte
tsn = dt/n;
t = linspace(ti(1),ti(2),n);

y(1) = y0;   % Startwert
option = optimset('TolX',1e-10);


for l=1:n-1
    % y(l+1) = y(l) + tsn * of(t(l+1),y(l+1));  impliziter Euler
    fh = @ (yr) (yr - y(l) - tsn * of(t(l+1),yr));
    y(l+1) = fzero(fh,y(l),option);
end