// resoud le pb dU/dt + dU^2/2/dx=0 (burgers)
// avec differents solveurs 

exec init.sci;
exec solveur.sci;

a=-1; b=1;
N=100;
dx=(b-a)/N;
x=[a:dx:b];
Tf=.8;

CFL=.8;



F=zeros(1,N+5);
div=zeros(1,N+5);

U=initialisation(a,b,N,x);
dt=CFL*dx/max(abs(U));

driver("X11")
xselect()
xbasc()
xset("pixmap",1)


t=0;
while t<Tf
 
 //conditions limites outflow
    U(1)=U(3);
    U(2)=U(3);
    U(N+4)=U(N+3);
    U(N+5)=U(N+3);

 
 //calcul des flux Fi+1/2

 for i=2:N+3,
   //F(i)=upwind(U,i);
   //F(i)=friedrich(U,i,dt,dx);
   F(i)=roe_O2(U,i);
   //F(i)=roe_llf(U,i,dt,dx);
 end

 for i=3:N+3,
   div(i)=dt/dx*(F(i)-F(i-1));
 end

   U1=U-div;

 //conditions limites outflow
    U1(1)=U1(3);
    U1(2)=U1(3);
    U1(N+4)=U1(N+3);
    U1(N+5)=U1(N+3);

 
 //calcul des flux Fi+1/2

 for i=2:N+3,
   //F(i)=upwind(U1,i);
   //F(i)=friedrich(U1,i,dt,dx);
   F(i)=roe_O2(U1,i);
   //F(i)=roe_llf(U1,i,dt,dx);
 end

 for i=3:N+3,
   div(i)=dt/dx*(F(i)-F(i-1));
 end

 U=0.5*(U+U1-div);





 t=t+dt;
 dt=CFL*dx/max(abs(U));

  xset("wwpc")
  plot2d(x,U(1,3:N+3));
  plot2d(x,U(1,3:N+3),-1);
  xset("wshow")
  xbasc()


end;






