function [t Y] = implicit_euler(F,Y0,T,h)
%backward euler method for systems/or scalar equation
%using built-in function fsolve
% input: 
%   F = @(t,X) [f_1(t,X(1),...,X(n)),...,f_n(t,X(1),...,X(n))]; function handle
%   Y0  row vector or a point, initial condition
%   T = [t0 tf]  row vector, t0 = starting time, tf = final time
%   h stepsize; optional argument, default h = 0.02
% 

%check input
if ~exist('h','var') || isempty(h)
    h = 0.02;
end

if (T(1) >= T(2))
    fprintf('Chyba: neplatny interval.\n');
    return
end

if h >= norm(T(2) - T(1),2)
    fprintf('Chyba: neplatne h.\n');
    return
end

%divide time interval
t = [T(1):h:T(2)]';
n = length(t);
if t(n) < T(2)
  t(n+1) = T(2);
elseif t(n) > T(2)
  t(n) = T(2);
end

%initialize for approximate values
N = length(t);
vars = length(Y0);
Y = zeros(N,vars);
Y(1,:) = Y0;

m = 0;


for i=1:(N-1)
    %first guess
    x0 = Y(i,:)';
    g = @(x) (x - h.*F(t(i+1),x) - Y(i,:)');    
    [extra x1] = evalc('fsolve(g,x0);');
    m = m + 1;
    Y(i+1,:) = x1';
end
 
printf('pocet reseni rovnice pro y_i+1 %f\n',m)
    
end
 

