gam := 1.4;
p[L] := 1;
u[L] := 0;
rho[L] := 1;
p[R] := 1;
u[R] := 1;
rho[R] := 1;
p1 = 1/2*(p[L]+p[R]);
A[L] := 2/((gam+1)*rho[L]);
A[R] := 2/((gam+1)*rho[R]);
B[L] := (gam-1)*p[L]/(gam+1);
B[R] := (gam-1)*p[R]/(gam+1);
a[L] := (gam*p[L]/rho[L])^(1/2);
a[R] := (gam*p[R]/rho[R])^(1/2);
if p1 > p[L] then
     f[L] := (p1-p[L])*(A[L]/(p1+B[L]))^(1/2);
     df[L] := (A[L]/(p1+B[L...