function [t,sol,x1,x2,x3,x4,err]=vari_lmm(l,h)
% Il programma calcola la soluzione approssimata dell'equazione differenziale
%   x'=l(x-sin(t))+cos(t)
%   x(0)=1
% utilizzando 4 metodi LMM. La soluzione esatta e' x(t)=exp(lt)+sin(t).
%
% La funzione restituisce 4 vettori:
% xi(j) contiene l'approssimazione della soluzione utilizzando il metodo i
% mentre err(j) contiene il max del valore assoluto dell'errore 
% globale di discretizzazione
%
% Con l=-10
% h=1.e-2, err(2:4) = 0.0402    0.0015    0.0011
%
% Con l=10, 
% h=1.e-2, err(2:4)= 1.0e+04 *[ 1.2925    0.0833    0.0689]
% h=1.e-4, err(2:4)=  218.8551, 0.0917, 0.0733
%
% Con l=-100
%  h=1.e-2, err(2:4)= 1.1353    0.1807    0.0479
% Osservare le forti oscillazioni dei grafici
% h=1.e-4, err(2:4)=3.7e-03  1.5e-05, 1.21e-05
%
% Proprieta' teoriche dei metodi:
% 1. Consistente, ma non 0-stabile (ha radici 1 e 3)
% 2. Consistente, ma "debolmente" 0-stabile (ha radici 1 e -1)
% 3. Esplicito. Consistente (di ordine 2) e 0-stabile ---> convergente di ordine 2
% 4. Implicito. Consistente (di ordine 2) e 0-stabile ---> convergente di ordine 2
x0=exp(l*h)+sin(h);
t=[0:h:1];
n=length(t);
x1=ones(1,n);
x2=ones(1,n);
x3=ones(1,n);
x4=ones(1,n);

% Il valore della soluzione nel secondo punto e' calcolato esattamente
x1(2)=x0;
x2(2)=x0;
x3(2)=x0;
x4(2)=x0;

    for i=2:n-1
        x1(i+1)=4*x1(i)-3*x1(i-1)-2*h*f(t(i-1),x1(i-1),l);
        x2(i+1)=x2(i-1)+2*h*f(t(i-1),x2(i-1),l);
        x3(i+1)=x3(i)+0.5*h*(3*f(t(i),x3(i),l)-f(t(i-1),x3(i-1),l));
        x4(i+1)=(1/(1-(2/3)*h*l))*((4/3)*x4(i)-(1/3)*x4(i-1)-(2/3)*h*l*sin(t(i+1))+(2/3)*h*cos(t(i+1))); 
    end


err=zeros(4,1);
sol=exp(l*t)+sin(t);
err(1)=max(abs(sol-x1));
   err(2)=max(abs(sol-x2));
   err(3)=max(abs(sol-x3));
   err(4)=max(abs(sol-x4));



function y=f(t,x,l)
y=l*(x-sin(t))+cos(t);
