program diffusion

      parameter( mx=500 )

      dimension u(mx+1,mx+1), r(mx+1,mx+1)

      pi = acos(-1.0)

c   *** Open files and write headers

      open(unit=1,file='residual.dat')
      open(unit=2,file='solution.dat')

      write(1,1)
1     format('#    nt       t       rms change')
      write(2,2)
2     format('#     x           u')

c   *** Prompt user for inputs

      print *, 'enter Nx, Ny, n'
      read(5,*) Nx, Ny, n

c   *** Set factors relating to the mesh

      Nxp1 = Nx+1
      Nyp1 = Ny+1

      dx = 1.0/float(Nx)
      dy = 1.0/float(Ny) 
      dx_inv  = 1.0/dx
      dx2_inv = 1.0/dx**2
      dy_inv  = 1.0/dy
      dy2_inv = 1.0/dy**2

      scale_xy = 1.0/(float(Nx-1)*float(Ny-1))

      dt = ?

c   *** Initialize the solution on the interior

      do j=2, Ny
         do i=2, Nx
            u(i,j) = 0.0
         end do
      end do

c   *** Initialize the solution on the boundaries

      do j=1, Ny+1
         u(1   ,j) = ?
         u(Nxp1,j) = ?
      end do
      arg = pi*(0.5+2.0*float(n))
      do i=1, Nx+1
         x = float(i-1)*dx
         u(i,1) = ?
         u(i,Nyp1) = ?
      end do

c   *** Initialize other parameters and enter the time stepping loop

      rms_old = 1.0e+10
      t = 0.0

      do nt=1, 10000000

c      *** Compute the right hand side

         do j=2, Ny
            do i=2, Nx
               d2udx2 = ?
               d2udy2 = ?
               r(i,j) = d2udx2 + d2udy2
            end do
         end do

c      *** Advance the solution one level in time.  Compute the rms change
c      *** to the solution at the same time.

         t = t + dt

         sum1 = 0.0
         do j=2, Ny
            do i=2, Nx
               du = dt*r(i,j)
               u(i,j) = u(i,j) + du
               sum1 = sum1 + du**2
            end do
         end do

         rms = sqrt( sum1*scale_xy )

c      *** Write the rms change only when it drops by an order of magnitude

         if( n .eq. 1 ) n_order = int(alog10(rms)) + 1

         m_order = int(alog10(rms))
         if( m_order .lt. n_order ) then
            n_order = m_order
            write(1,10) nt, t, rms
10          format(i8, 1p,2e12.4)
         end if

c      *** Check for convergence

         if( rms .lt. 1.0e-7 ) goto 50

c      *** Make sure the solution is not diverging

         if( rms .gt. rms_old ) then
            print *, 'solution diverging'
            stop
         end if

         rms_old = rms

      end do

c   *** If flow gets here, the solution was not converged

      print *, 'WARNING: the convergence condition was not met'

50    continue

c   *** Write the final solution

      j = ?
      do i=1, Nx+1
         x = float(i-1)*dx
         write(2,20) x, u(i,j)
20       format(1p,2e12.4)
      end do

      stop
      end