# codigo para resolver q_t + ( umax q (1-q) )_x = 0
#dimensionamentos
n = 200;
qa = zeros(n,1);
q = qa;
f = zeros(n+1,1);
x = q;
#dados
t = 0;
dt = 0.001;
nt = 600;
dx = 1./n;
qezq = 0.3;                              
qder = 0;

#condicao inicial
for i=1:n
x(i)=(i-0.5)*dx;
q(i)=0.3;
endfor
plot(x,q,"-"); pause

#velocidades maximas
for i=1:n+1
 xx=(i-1)*dx;
 umax(i)=1-0.5*exp(-(xx-0.5)^2/(0.05)^2);
 dumax(i)=-0.5*exp(-(xx-0.5)^2/(0.05)^2)*(-2)*(xx-0.5)/(0.05)^2;
endfor

#codigo
for m=1:nt
# atualizamos
        qa = q;
# calculo dos fluxos
for i=2:n
  xx = 0.5*(x(i-1)+x(i));
  qm = (qa(i-1)+qa(i))/2; #media
  beta = (1-2*qm)*umax(i);
#  if (beta>0) qq = qa(i-1); else qq=qa(i); endif 
                 #full upwind, Godunov 
  fr=umax(i)*q(i)*(1-q(i));
  fl=umax(i)*q(i-1)*(1-q(i-1));
  qq=qm-dt/2/dx*(fr-fl);
                 #richtmyer-lax-wendroff 
  f(i)=umax(i)*(1-qq)*qq; 
endfor

#ezquerda
qm=(qezq+qa(1))/2;
beta = (1-2*qm)*umax(1);
if (beta>0) qq = qezq; else qq=qa(1); endif 
                 #full upwind, Godunov 
f(1)=umax(1)*(1-qq)*qq; 

#direita
qm=(qder+qa(n))/2;
beta = (1-2*qm)*umax(n+1);
if (beta>0) qq = qa(n); else qq=qder; endif 
                 #full upwind, Godunov 
f(n+1)=umax(n+1)*(1-qq)*qq; 

# calculo do q
for i=1:n
	q(i)=qa(i)-dt/dx*(f(i+1)-f(i));
endfor
# graficamos
	plot(x,q,"-+"); pause
 t=t+dt;
endfor