      program heatMOL
      INTEGER NSTEP,NVAR
      PARAMETER(NVAR=100,NSTEP=3001)
      INTEGER i,j,N
      REAL xx(NSTEP),x1,x2,y(NVAR,NSTEP),v(NVAR)
      REAL sigma,b,h,k,r,T,L,M
      COMMON /parameters/ sigma,b,N,r
      COMMON /path/ xx,y
      EXTERNAL derivs
 

      sigma=10.0
      b=8.0/3.0
      
      L=10.0
      T=10.0
      N=NVAR
      M=3000
      h=L/( N+1.0);
      k=T/ M;
      r=1.0/(h*h); 
      
      x1=0.0


      do i=1,N
         v(i)=100
      enddo
      v(1)=0
      v(N)=0

      
      call rkdumb(v,N,x1,T,NSTEP-1,derivs)

      open (unit=10,file="tempMOL.dat",status="unknown")
      do 1 i=1,N
         x1=i*h
      write(10,*) x1, y(i,NSTEP-1)
1     continue

      close(10)

      stop
      end

      subroutine derivs(x,y,dy)

      real x,y(*),dy(*),sigma,r,b
      integer i
      common /parameters/ sigma,b,N,r
c     y(0)=0
      dy(1)=r*(-2*y(1)+y(2))
      do i=2,N-1
      dy(i)=r*(y(i-1)-2*y(i)+y(i+1))
      enddo
      dy(N)=r*(y(N-1)-2*y(N))
c     y(N+1)=0
      return
      end


      SUBROUTINE rkdumb(vstart,nvar,x1,x2,nstep,derivs)
      INTEGER nstep,nvar,NMAX,NSTPMX
      PARAMETER (NMAX=100,NSTPMX=3001)
      REAL x1,x2,vstart(nvar),xx(NSTPMX),y(NMAX,NSTPMX)
      EXTERNAL derivs
      COMMON /path/ xx,y
CU    USES rk4
      INTEGER i,k
      REAL h,x,dv(NMAX),v(NMAX)
      do 11 i=1,nvar
        v(i)=vstart(i)
        y(i,1)=v(i)
11    continue
      xx(1)=x1
      x=x1
      h=(x2-x1)/nstep
      do 13 k=1,nstep
        call derivs(x,v,dv)
        call rk4(v,dv,nvar,x,h,v,derivs)
        if(x+h.eq.x)pause 'stepsize not significant in rkdumb'
        x=x+h
        xx(k+1)=x
        do 12 i=1,nvar
          y(i,k+1)=v(i)
12      continue
13    continue
      return
      END
