c*************************************************************
c*
c*  Implicit method for solution of heat conduction problem...
c*
c*      u(xx) = u(t) + f(x,t)  for 0<x<xR  and 0<t<tmax
c*
c*      where u(x,0)=g(x), u(0,t)=uL, and u(L,t)=uR
c*
c*************************************************************
      implicit real*8(a-h,o-z)
      real*8 lam, L
      parameter(L=1.0,uL=0.0,uR=0.0,nx=40,nt=50,tmax=0.5)
      parameter(a=0.2,b=0.6)
      dimension x(0:nx+1), t(0:nt), u(0:nx+1)
      dimension w(nx),v(nx),z(nx)
      f(x,t)=0.0
c      g(x)=x*(1.0-x)
c      g(x)=1.0
      g(x)=step(x-a)-step(x-b)
      dx=L/(nx+1)
      dt=tmax/nt
      lam=dt/(dx*dx)	
      xa=-1.0-2.0*lam
      t(0)=0.0
      open(unit=6,file="data")
      write(6,103)
  103 format('U:=[')
      do 2 ix=0,nx+1
        x(ix)=dx*ix
    2   u(ix)=g(x(ix))
      write(6,104)
      write(6,100)  (t(0),x(ix),u(ix),ix=0,nx)
	write(6,107)  t(0),x(nx+1),u(nx+1)
      write(6,105)
	u(0)=uL
	u(nx+1)=uR
      do 40 it=1,nt
        t(it)=dt*it
	  w(1) = xa
        v(1) = lam/w(1)
        rr=-u(1)+dt*f(x(1),t(it))
	  z(1)=rr/w(1)
        do 20  jj = 2, nx
	   w(jj)=xa-lam*v(jj-1)
	   v(jj) = lam/w(jj)
           rr=-u(jj)+dt*f(x(jj),t(ij))
   20      z(jj) = (rr-lam*z(jj-1))/w(jj)
	u(nx) = z(nx)
	do 30  jj = nx - 1, 1, -1
   30      u(jj) = z(jj) -v(jj)*u(jj+1)
      write(6,104)
  104 format('[')
      write(6,100)  (t(it),x(ix),u(ix),ix=0,nx)
  100 format('[',f8.4,',',f8.4,',',e12.4,'],')
      write(6,107)  t(it),x(nx+1),u(nx+1)
  107 format('[',f8.4,',',f8.4,',',e12.4,']')
      if(it.ne.nt) write(6,105)
  105 format('],')
   40 continue
      write(6,106)
  106 format(']];')
      close(unit=6)
      end

      function step(x)
      implicit real*8(a-h,o-z)
      step=0.0
      if(x.ge.0.0) step=1.0
      return
      end