function [t Y] = implicit_euler_newton(F,Y0,T,h,dF)
%backward euler method for systems/or scalar equation
%using Newton's method
% 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; default h = 0.02
%   dF optional argument, funtion handle, Jacobi of F


%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;



for i=1:(N-1)
    %g(y_{i+1}) = y_{i+1} - hF(t_{i+1},y_{i+1}) - Y_i
    g = @(t,x) (x - h.*F(t,x)- Y(i,:)');
    %newton method set up
    tol = 1e-3;
    max_iters = 50;
    c = true;
    k = 1;
    
    %initial guess
    x0 = Y(i,:)';
    t0 = t(i);
    %check if dF provided
    if ~exist('dF','var')
      J = jacobian(F,t0,x0,h);
    else
      J = eye(length(x0)) - h.*dF(t0,x0);
    end
    %no inverse version
    delta = J \ (-g(t0,x0));
    x1 = x0 + delta;
    %x1 = x0 - inv(eye(length(x0)) - h.*dF(t0,x0))*(g(t0,x0));
    %x1 = x0 - (inv(eye(length(x0)) - h.*dF(t0,x0))*(x0 - h.*F(t0,x0) - Y(i,:)'))
    
    while norm(x1-x0,inf) > tol && c
      x0 = x1;
      %check if dF provided
      if ~exist('dF','var')
        J = jacobian(F,t0,x0,h);
      else
        J = eye(length(x0)) - h.*dF(t0,x0);
      end
            
      delta = J \ (-g(t0,x0));
      x1 = x0 + delta;
      %x1 = x0 - (inv(eye(length(x0)) - h.*dF(t0,x0))*(x0 - h.*F(t0,x0) - Y(i,:)'))
      %x1 = x0 - inv(eye(length(x0)) - h.*dF(t0,x0))*(g(t0,x0));
      
      k = k + 1;
      if k == max_iters
        c = false;
        printf('newton reached max_iters\n')
      end
     end
     Y(i+1,:) = x1';      
end   
  
    
end
 
