function [ tout, yout ] = impEuler(f, tspan, y0, h, tol)
% Implementierung des impliziten Eulerverfahrens mit
% Fixpunktiteration fuer die Loesung des Gleichungssystems

%MATLABALL:NO_EXEC
    
	if nargin < 1
		% Keine Argumente -> automatisierter Aufruf, beenden
		disp('Keine Argumente uebergeben, beende.');
	end
	
	if nargin < 5
		tol = h * 1e-2;
	end
	
	N = round(diff(tspan) / h);
	h = diff(tspan) / N;
	
	tout = linspace(tspan(1), tspan(2), N + 1);
	yout = zeros(length(y0),            N + 1);
	yout(:, 1) = y0;
	
	t = tspan(1);
	y = y0;
	
	for kk = 1:N
		% Predict with an explicit Euler method
		Y1 = y + h * f(t, y);
		
		Phi = @(ynp) y + h * f(t + h, ynp);
		
		PhiY1 = Phi(Y1);
		
		while (norm(PhiY1 - Y1) > tol)
			Y1    = PhiY1;
			PhiY1 = Phi(Y1);
		end
		
		y = y + h * f(t + h, Y1);
		t = t + h;
		
		yout(:, kk + 1) = y;
	end
end
