function [t,y] = euler_back_secants(t0,y0,t_end,h,fcn,maxit,tol)

%----------------------------------------------------
% VARIANT OF A MATLAB CODE TAKEN FROM 
% K. Atkinson, W. Han e D. Steward, 
% Numerical Solution of Ordinary Differential Equations, 
% programs for the book, John Wiley and Sons, 2009, 
% http://www.math.uiowa.edu/NumericalAnalysisODE/
%----------------------------------------------------
% Solve the initial value problem
%    y' = f(t,y),  t0 <= t <= b,  y(t0)=y0
% Use the backward Euler method with a stepsize of h.  The user must
% supply an m-file to define the derivative f, with some name,
% say 'deriv.m', and a first line of the form
%   function ans=deriv(t,y)
% tol is the user supplied bound on the difference between successive
% values of the backward Euler iteration.
% Maxit is the maximum number of iterations performed
% by the secants method.
% A sample call would be
%   [t,z]=euler_back_secants(t0,z0,b,delta,'deriv',maxit,1e-6)
%----------------------------------------------------
% Output:
% The routine euler_back will return two vectors, t and y.
% The vector t will contain the node points
%   t(1)=t0, t(j)=t0+(j-1)*h, j=1,2,...,N
% with
%   t(N) <= t_end,  t(N)+h > t_end
% The vector y will contain the estimates of the solution Y
% at the node points in t.
%----------------------------------------------------

% Initialize.
n = fix((t_end-t0)/h) + 1;
t = linspace(t0,t0+(n-1)*h,n)';
y = zeros(n,1);
y(1) = y0;
% advancing
for i=2:n
    % forward Euler estimate
    yt1 = y(i-1) + h*feval(fcn,t(i-1),y(i-1));
    
    % fixed point iteration
    yt2 = y(i-1) + h*feval(fcn,t(i),yt1);
    
    % secant iteration
    k=1;
    yy(1)=min(yt1,yt2); yy(2)=max(yt1,yt2);
    step_err(1)=yy(2)-yy(1);

    while (abs(step_err(k)) > tol) & (k < maxit)

        k=k+1;

        fx_old=yy(k-1)-y(i-1)-h*feval(fcn,t(i),yy(k-1));
        fx_new=yy(k)-y(i-1)-h*feval(fcn,t(i),yy(k));
        dfx=(fx_new-fx_old)/(yy(k)-yy(k-1));

        step_err_k=-fx_new/dfx;
        step_err=[step_err; step_err_k];

        yy(k+1)=yy(k)+step_err(end);
        
    end
    
    step_err=[];
    
    if k == maxit 
    fprintf('\n \t Results may be not correct.');
    end

    y(i) = yy(end);
end