      subroutine intpar(neqn,nout,nsteps,t0,tf,abserr,relerr)
C
C Subroutine intpar sets the parameters to control the
C integration of the linear PDE problem
C
C Double precision coding is used 
      implicit double precision(a-h,o-z)
C
C Number of ODEs
      neqn=21
C
C Number of output points
      nout=6
C
C Maximum number of steps in the interval t0 to tf
      nsteps=250
C
C Initial, final values of the independent variable
      t0=0.0d0
      tf=0.2d0
C
C Error tolerances
      abserr=1.0d-05
      relerr=1.0d-05
      return
C
C End of intpar
      end
      subroutine par(xl,xu,pi)
C
C Subroutine par sets the parameters for the linear
C PDE problem
C
C Double precision coding is used
      implicit double precision(a-h,o-z)
C
C Problem parameters
      xl=0.0d0
      xu=1.0d0
      pi=4.0d0*datan(1.0d0)
      return
C
C End of par
      end
      subroutine inital(neqn,t,u0)
C
C Subroutine inital sets the initial condition vector 
C for the linear PDE problem
C
C Double precision coding is used
      implicit double precision(a-h,o-z)
C
C Size the arrays
      dimension u0(neqn)
C
C Problem parameters
      call par(xl,xu,pi)
C
C Initial condition
      do i=1,neqn
        x=xl+dfloat(i-1)/dfloat(neqn-1)*(xu-xl)
        u0(i)=dsin(pi*x)
      end do
      return
C
C End of inital
      end
      subroutine derv(neqn,t,u,ut)
C
C Subroutine derv computes the derivative vector
C of the linear PDE problem
C
C Double precision coding is used
      implicit double precision(a-h,l,o-z)
C
C Size the arrays
      dimension u(neqn), ut(neqn)
C
C Problem parameters
      call par(xl,xu,pi)
C
C BC at x = 0
      ut(1)=0.0d0
C
C BC at x = 1
      ut(neqn)=0.0d0
C
C Interior points
      dx=(xu-xl)/dfloat(neqn-1)
      dxs=dx*dx
      do i=2,neqn-1
        ut(i)=(u(i+1)-2.0d0*u(i)+u(i-1))/dxs
      end do
      return
C
C End of derv
      end
      subroutine fprint(no,ncase,neqn,t,u)
C
C Subroutine fprint displays the numerical and
C exact solutions to the linear PDE problem
C
C Double precision coding is used
      implicit double precision(a-h,l,o-z)
C
C Size the arrays
      dimension u(neqn)
C
C Problem parameters
      call par(xl,xu,pi)
C
C Print a heading for the solution at t = 0
      if(t.le.0.0d0)then
C
C       Label for ODE integrator
C
C         Fixed step modfied Euler
          if(ncase.eq.1)then
            write(no,11)
 11         format(//,' euler2a integrator')
C
C         Variable step modified Euler
          else if(ncase.eq.2)then
            write(no,12)
 12         format(//,' euler2b integrator')
C
C         Fixed step classical fourth order RK
          else if(ncase.eq.3)then
            write(no,13)
 13         format(//,' rkc4a integrator')
C
C         Variable step classical fourth order RK
          else if(ncase.eq.4)then
            write(no,14)
 14         format(//,' rkc4b integrator')
C
C         Fixed step RK Fehlberg 45
          else if(ncase.eq.5)then
            write(no,15)
 15         format(//,' rkf45a integrator')
C
C         Variable step RK Fehlberg 45
          else if(ncase.eq.6)then
            write(no,16)
 16         format(//,' rkf45b integrator')
          end if
C
C Heading
      write(no,4)ncase,neqn
 4    format(/,' ncase = ',i1,'  neqn = ',i2)
      write(no,2)
 2    format(/,4x,'t',4x,'u1(num)',5x,'u1(ex)',8x,'diff1',/)
C
C End of t = 0 heading
      end if
C
C Numerical and analytical solution output 
C
C     Midpoint value of x 
      x=(xu-xl)/2.0d0  
C
C     Analytical solution at midpoint 	
      ue=dexp(-pi*pi*t)*dsin(pi*x)
C  	
C     Grid index of midpoint 
      im=dfloat((neqn+1)/2)
C	
C     Display the numerical and exact solutions, and their difference 
      write(no,3)t,u(im),ue,u(im)-ue
 3    format(f5.2,2f11.6,e13.4)
      return
C
C End of fprint
      end
