#dimensioning
n = 600;
qa = ones(2,n);
q = qa;
f = zeros(2,n+1);
x = zeros(1,n);
#dados
u0 = 0.;
h0 = 1.;
g=1.;
c0=sqrt(g*h0);
cfl=0.9;
nt = ceil(n/cfl);
t = 0;
dt = 1./nt;
dx = 1./n;
for i=1:n;
 x(i)=(i-0.5)*dx;
 eta=2.*exp(-(x(i)-0.4)^2/0.002);
 q(1,i)=eta;
 q(2,i)=0;
endfor
#code----------------------------------
plot(x,q(1,:),"-"); axis([0 1 -1 2]); pause
A=[0 1;-u0*u0+g*h0 2*u0];
R=[1 1;u0-c0 u0+c0];
D=[u0-c0 0;0 u0+c0];
Dmenos=min(D,0);
Imenos=abs(sign(Dmenos));
Mmenos=R*Imenos*inv(R);
for m=1:nt
# update
        qa = q;
# fluxes
i=2:n;
qint_godunov=qa(:,i-1)+Mmenos*(qa(:,i)-qa(:,i-1));
f(:,i)=(A*qint_godunov);
#boundary conditions
w11=1/(2*c0)*((u0+c0)*qa(1,1)-qa(2,1));
#closed left
hezq=2*c0*w11/(u0+c0);
qezq=[hezq; 0];
f(:,1)=A*qezq;
#closed right
w2n=1/(2*c0)*(-(u0-c0)*qa(1,n)+qa(2,n));
hder=2*c0*w2n/(c0-u0);
qder=[hder;0];
f(:,n+1)=A*qder;

# computing q
i=1:n;
q(:,i)=qa(:,i)-dt/dx*(f(:,i+1)-f(:,i));
# plot
if ((mod(m,30)==0))
plot(x,q(1,:),"-+"); axis([0 1 -1 2]); pause
endif
t=t+dt;
endfor