C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C 1-dimensional Diffusion Equation with simple
C Dirichlet boundary conditions u(0,t)=u(1,t)=0
C 0<= x <= 1 and 0<= t <= tf
C 
C We set initial condition u(x,t=0) that satisfies 
C the given boundary conditions. 
C Nx is the number of points in spatial lattice:
C x = 0 + (j-1)*dx, j=1,...,Nx and dx = (1-0)/(Nx-1)
C Nt is the number of points in temporal lattice:
C t = 0 + (j-1)*dt, j=1,...,Nt and dt = (tf-0)/(Nt-1)
C
C u(x,0) = sin(pi*x) tested against analytical solution
C u(x,t) = sin(pi*x)*exp(-pi*pi*t)
C
C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

      program diffusion_1d
      implicit none
      integer P            ! Max no of points
      real*8  PI
      parameter (P=100000,PI=3.1415926535897932D0)
      real*8 u(P), d2udx2(P)
      real*8 t,x,dx,dt,tf,courant
      integer Nx,Nt,i,j
C --- Input:
      print *, '# Enter: Nx, Nt, tf: (P= ',P,' Nx must be < P)'
      read(5,*) Nx,Nt,tf
      if(Nx .ge. P) stop 'Nx >= P'
      if(Nx .le. 3) stop 'Nx <= 3'
      if(Nt .le. 2) stop 'Nt <= 2'
C --- Initialize:
      dx      = 1.0D0/(Nx-1)
      dt      = tf   /(Nt-1)
      courant = dt/dx**2
      print * ,'# 1d Diffusion Equation: 0<=x<=1, 0<=t<=tf'
      print * ,'# dx= ',dx,' dt= ',dt,' tf= ', tf
      print * ,'# Nx= ',Nx,' Nt= ',Nt
      print * ,'# Courant Number= ',courant
      if(courant .gt. 0.5D0) print *,'# WARNING: courant > 0.5'
      open(unit=11,file='d.dat') ! data file
C --- Initial condition at t=0 ------------------------------
C     u(x,0) = sin( pi x) e^{-pi^2 t}
      do i= 1, Nx
       x    = (i-1)*dx
       u(i) = sin(PI*x)
c       u(i) = sin(PI*x*2.0D0)
      enddo
      u(1)  = 0.0d0
      u(Nx) = 0.0d0
      do i= 1,Nx
       x    = (i-1)*dx
       write(11,*) 0.0D0, x, u(i)
      enddo
      write(11,*)' '

C ----------------------------------------------------------

C --- Calculate time evolution:
      do j=2,Nt
       t = (j-1)*dt
C ----- second derivative:
       do i=2,Nx-1
        d2udx2(i) = courant*(u(i+1)-2.0D0*u(i)+u(i-1))
       enddo
C ----- update:
       do i=2,Nx-1
        u(i) = u(i) + d2udx2(i)
       enddo
       do i=1,Nx
        x =  (i-1)*dx
        write(11,*) t, x, u(i)
       enddo
       write(11,*)' '
       
      enddo ! do j=2,Nt

      close(11)
      end
