// resoud le pb dU/dt +c dU/dx=0
// avec le schema upwind ordre 1

function [U]=initialisation(a,b,N,c,x)
   U=zeros(1,N+5);

  //une gaussienne centree en 0
   y=exp(-x.*x*5);
   U(1,3:N+3)=y;

 //une fonction discontinue en 0
   //y=sign(x);
   //U(1,3:N+3)=y;

endfunction;


a=-5; b=5;
N=200;
dx=(b-a)/N;
x=[a:dx:b];
Tf=90.;

CFL=.9;

c=0.3;
cplus=max(c,0);
cmoins=min(c,0);

dt=CFL*dx/abs(c);

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

U=initialisation(a,b,N,c);



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

 
 //calcul des flux Fi+1/2

 for i=2:N+3,
   F(i)=cplus*U(i)+cmoins*U(i+1);
 end

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

 U=U-div;
 t=t+dt;

xset("window",0)
  drawlater()
  clf();

  plot2d(x,U(1,3:N+3));
  
 drawnow()
end;






