function [t,y,flag]=eulero_implicito_heun(t0,y0,t_end,h,fcn,tol)

% Metodo di Eulero implicito

flag=0; max_count=100;

n = fix((t_end-t0)/h) + 1;
t = linspace(t0,t0+(n-1)*h,n)';
y = zeros(n,1); y(1) = y0;

for i=2:n
    % valutazione soluzione in t(i)
    f1=feval(fcn,t(i-1),y(i-1)); % EE
    f2=feval(fcn,t(i),y(i-1)+h*f1); % HEUN
    yt1 = y(i-1) + (h/2)*(f1+f2);
    
    count = 0;  diff = 1;
    while diff(end) > tol && count < max_count % iterazioni punto fisso
        yt2 = y(i-1) + h*feval(fcn,t(i),yt1);
        diff(count+1) = abs(yt2-yt1);
        yt1 = yt2;
        count = count +1;
    end
    
    y(i) = yt2;
    
    if count >= max_count, flag=1; end
    
end