
function esercizio3

t0=0;
y0=0;
t_end=10;
fcn=@(x,y) 0*x+(cos(y)).^2;
fsol=@(x) atan(x);

% nb: mettiamo piu' punti rispetto all'esercizio
NV=[8 16 32 64 128 256];

for k=1:length(NV)

    N=NV(k);

    h=(t_end-t0)/(N-1);

    % soluzione numerica (t,u: vettori colonna)
    % [t,u]=adams_bashforth_3passi(t0,y0,t_end,h,fcn,fsol);
    [t,u]=adams_bashforth_4passi(t0,y0,t_end,h,fcn,fsol);

    % soluzione esatta
    y=fsol(t');

    % errore assoluto
    errV=abs(u-y);
    abserr(k)=max(errV);

    % plot k-simo grafico
    figure(k);
    semilogy(t,errV);
    stringa_titolo=strcat('N: ',num2str(N)); 
    title(stringa_titolo);
 
end

figure(length(NV)+1);
semilogy(NV,abserr,'r-o','LineWidth',2,...
                       'MarkerEdgeColor','k',...
                       'MarkerFaceColor','g',...
                       'MarkerSize',10);
title('Norma infinito errori');




function [t,y]=adams_bashforth_3passi(t0,y0,t_end,h,fcn,fsol)

% Oggetto: metodo di Adams-Bashforth a 3 passi.
% t,y: vettori colonna di ascisse e soluzione numerica.

% ascisse
t=t0:h:t_end;

% passi iniziali
t_end_start=t(3);

% [t_start,y_start]=eulero_esplicito(t0,y0,t_end_start,h,fcn);
% [t_start,y_start,flag]=eulero_implicito(t0,y0,t_end_start,h,fcn,10^(-12));
% [t_start,y_start,flag]=crank_nicolson(t0,y0,t_end_start,h,fcn,10^(-12));
y_start=feval(fsol,(t0:h:t_end_start)');

y=y_start;

% passi successivi
N=length(t);
for k=4:N
    fnm3=feval(fcn,t(k-3),y(k-3));
    fnm2=feval(fcn,t(k-2),y(k-2));
    fnm1=feval(fcn,t(k-1),y(k-1));
    y(k)=y(k-1)+(h/12)*(23*fnm1-16*fnm2+5*fnm3);
end






function [t,y]=adams_bashforth_4passi(t0,y0,t_end,h,fcn,fsol)

% Oggetto: metodo di Adams-Bashforth a 4 passi.
% t,y: vettori colonna di ascisse e soluzione numerica.

% ascisse
t=t0:h:t_end;

% passi iniziali
t_end_start=t(4);
% [t_start,y_start]=eulero_esplicito(t0,y0,t_end_start,h,fcn);
% [t_start,y_start,flag]=eulero_implicito(t0,y0,t_end_start,h,fcn,10^(-12));
[t_start,y_start,flag]=crank_nicolson(t0,y0,t_end_start,h,fcn,10^(-12));
% y_start=feval(fsol,t0:h:t_end_start);

y=y_start;

% passi successivi
N=length(t);
for k=5:N
    fnm4=feval(fcn,t(k-4),y(k-4));
    fnm3=feval(fcn,t(k-3),y(k-3));
    fnm2=feval(fcn,t(k-2),y(k-2));
    fnm1=feval(fcn,t(k-1),y(k-1));
    y(k)=y(k-1)+(h/24)*(55*fnm1-59*fnm2+37*fnm3-9*fnm4);
end


