function [xx,tt,u]=hyper(xspan,tspan,u0,ul,...
                          scheme,cfl,deltax,deltat)
% HYPER approssimazione eq scalare iperbolica, a>0
% [XX,TT,U]=HYPER(XSPAN,TSPAN,U0,UL,SCHEME,CFL,...
%                 DELTAX,DELTAT)
% risolve l'equazione differenziale iperbolica scalare
%       DU/DT+ A * DU/DX=0
% in (XSPAN(1),XSPAN(2))x(TSPAN(1),TSPAN(2))
% con condizione iniziale  U(X,0)=U0(X) e
% condizione al bordo U(T)=UL(T) assegnata in XSPAN(1)
% con vari schemi alle differenze finite.
% scheme = 1 Lax - Friedrichs
%          2 Lax - Wendroff
%          3 Upwind
% La velocita' di propagazione A non e' richiesta
% esplicitamente, essendo CFL = A * DELTAT / DELTAX
% In output XX e' il vettore della discretizzazione
% in x; TT e' il vettore della discretizzazione in t
% U e' una matrice che contiene la soluzione numerica:
% U(n,:) contiene la sol all'istante temporale TT(n)
% U0 e UL possono essere inline o anonymous function o
% function definite tramite M-file.

Nt=(tspan(2)-tspan(1))/deltat+1;
tt=linspace(tspan(1),tspan(2),Nt);
Nx=(xspan(2)-xspan(1))/deltax+1;
xx=linspace(xspan(1),xspan(2),Nx);
u=zeros(Nt,Nx); cfl2=cfl*0.5; cfl21=1-cfl^2;
cflp1=cfl+1; cflm1=cfl-1;
u(1,:)=feval(u0,xx);
for n=1:Nt-1
 u(n+1,1)=feval(ul,tt(n+1));
 if scheme == 1
% Lax Friedrichs
    for j=2:Nx-1
      u(n+1,j)=0.5*(-cflm1*u(n,j+1)+cflp1*u(n,j-1));
    end
    j=Nx;
    u(n+1,j)=0.5*(-cflm1*(2*u(n,j)-u(n,j-1))+...
        cflp1*u(n,j-1));
 elseif scheme == 2
% Lax Wendroff
    for j=2:Nx-1
     u(n+1,j)=cfl21*u(n,j)+...
         cfl2*(cflm1*u(n,j+1)+cflp1*u(n,j-1));
    end
    j=Nx;
    u(n+1,j)=cfl21*u(n,j)+...
     cfl2*(cflm1*(2*u(n,j)-u(n,j-1))+cflp1*u(n,j-1));
 elseif scheme ==3
% Upwind
   for j=2:Nx
        u(n+1,j)=-cflm1*u(n,j)+cfl*u(n,j-1);
    end
 end
end