c code to solve, and then plot, the solution of the nonlinear BVP
c        u_t + u*u_x = ep*u_xx     for xL < x < xr  
c where 
c          u(xl,t) = a0 ,  u(xR,t) = b0,  u(x,0) = phi(x)

c    gfortran burgers.f	./a.out
c    gfortran burgers.f -O3 -ftree-vectorize -ffast-math

	implicit real*8(a-h,o-z)
	parameter( xL = 0.0 , a0 = 1.0 , xR = 1.0 , b0 = -1.0 )
	parameter( nx = 10000 , nt = 10000, tmax = 10000, ep = 0.01 )
	dimension  y(0:nx+1), a(nx+1), c(nx+1), u(nx+1), v(nx+1)
	dimension  x(0:nx+1), uu(0:nx+1), q(0:nx+1)
	f(xf,yf,zf, qf, ep, dt) = (yf*zf+yf/dt)/ep+qf
	fy(xf,yf,zf, qf, ep, dt) = (zf+1/dt)/ep
	fz(xf,yf,zf, qf, ep, dt) = yf/ep

	error = 0.000001

c  x(0)=xL  x(nx+1)=xR
	dx=(xR-xL)/(nx+1)
	dxx=dx*dx
	do 2 ix=0,nx+1
	  x(ix)=xL+dx*ix
	  uu(ix)=phi(x(ix))
    2 	continue

c t(1)=0  t(nt)=tmax
	dt=tmax/(nt-1)
	do 80 it=2,nt
	  t=dt*(it-1)

	do 10 ix=0,nx+1
	  q(ix)=-uu(ix)/(dt*ep)
	  y(ix)=uu(ix)
   10	continue

	err=1

	do while (err > error)

	do 20 j = 1, nx
	  z = (y(j+1) - y(j-1))/(2.0*dx)
	  a(j) = 2.0 + dxx*fy(x(j), y(j), z, q(j), ep, dt)
	  c(j) = -1.0 - 0.5*dx*fz(x(j), y(j), z, q(j), ep, dt)
	  v(j) = -(2.0*y(j)-y(j+1)-y(j-1)+dxx*f(x(j),y(j),z,q(j),ep,dt))
   20	continue
	v(1) = v(1)/a(1)
	u(1) = - (2.0 + c(1))/a(1)
	do 30 j = 2, nx
	  xxl = a(j) - c(j)*u(j-1)
	  v(j) = (v(j) - c(j)*v(j-1))/xxl
	  u(j) = - (2.0 + c(j))/xxl
   30	continue
	vv = v(nx)
	y(nx) = y(nx) + vv
	err = abs(vv)
	do 40 jj = nx - 1, 1, -1
	  vv = v(jj) - u(jj)*vv
	  err = max(err, abs(vv))
	  y(jj) = y(jj) + vv
   40	continue

	enddo


	do 60 ix=0,nx+1
	  uu(ix)=y(ix)
   60	continue

   80	continue

	open(unit=6,file="data_10000.txt")
c	  write(6,604)  tmax
	do 82 ix=0,nx+1
	  write(6,604)  x(ix), y(ix)
  604 format(e14.5,e14.5)
   82	continue
	close(unit=6)

	stop
	end

	function phi(x)
	implicit real*8(a-h,o-z)
	x0=0.25
	a=1
	if (x