function [t,x]=df_es2(L,N,x0,itmax,tol)
% Metodo alle differenze finite per il problema
% Exercise 7.4, page 187 Ascher-Petzold
% Si applica Newton
% Shooting NON funziona per questo esempio

% scegli ad esempio N=100; x0=ones(2*(N+1),1);
% [t,x]=df_es2(100,N,x0,15,1.e-10);
t=linspace(0,L,N+1);
h=1/N;
m=2;
S=zeros(m,m,N);
R=zeros(m,m,N);
f=zeros(m*(N+1),1);
x=zeros(m*(N+1),1);
B0=zeros(2);
B1=zeros(2);
B0(1,2)=1;
B1(2,1)=1;
H=spalloc(m*(N+1),m*(N+1),2*m^2*(N+1));
H(N*m+1:(N+1)*m,1:m)=B0;
H(N*m+1:(N+1)*m,N*m+1:(N+1)*m)=B1;
y=x0;
w=y;

for k=1:itmax
% un passo di Newton: x=y-inv(H)*w

% Aggiorna la matrice H e il vettore w
for n=1:N
    t1=t(n)+h/2;
    v1=(y(2*n-1)+y(2*n+1))/2;
    v2=(y(2*n)+y(2*n+2))/2;
    S(:,:,n)=[-1, -h/2;
        h*(t1*v1-1/2), -1+2*h/t1];
    
    R(:,:,n)=[1, -h/2;
        h*(t1*v1-1/2), 1+2*h/t1];
    
    w(2*n-1)=y(2*n+1)-y(2*n-1)-h*v2;
    
    w(2*n)=y(2*n+2)-y(2*n)+h*(4/t1*v2+(t1*v1-1)*v1);
end
w(2*N+1)=y(2);
w(2*N+2)=y(2*N+1);

for n=1:N
    H((n-1)*m+1:n*m,(n-1)*m+1:n*m)=S(:,:,n);
    H((n-1)*m+1:n*m,n*m+1:(n+1)*m)=R(:,:,n);
end

% Correzione di Newton
inc=H\w;

noinc=max(abs(inc));
fprintf('Iterazione = %4.0f, Incremento di Newton = %1.4e \n', k, noinc)
% passo di Newton
x=y-inc;
    if(noinc<tol)break
    end
y=x;
end
x=x(1:2:end);


    