      program pdenonm
C
C Main program pdenonm computes the numerical 
C solution to a nonlinear PDE by six integrators
C
C     Double precision coding is used
      implicit double precision(a-h,o-z)
C
C Size arrays
      parameter(neq=500)
      dimension u0(neq),u(neq)
C
C Open a file for output
      no=2
      open(no,file='pdenonfor.out')
C
C Step through six integrators
      do ncase=1,6
C
C       Integration parameters
        call intpar(neqn,nout,nsteps,t0,tf,abserr,relerr)
C
C       Initial condition 
        call inital(neqn,t0,u0)
C
C       Output time
        tp=tf-t0
C
C       Step through nout grid points
        do j=1,nout
C
C         Print solution
          call fprint(no,ncase,neqn,t0,u0)
C
C         Select ODE integrator
C
C         Fixed step modified Euler integrator
          if(ncase.eq.1)then
            call euler2a(neqn,t0,tf,u0,nsteps,u)
          end if
C
C         Variable step modified Euler integrator
          if(ncase.eq.2)then
            call euler2b(neqn,t0,tf,u0,nsteps,abserr,relerr,u)
          end if
C
C         Fixed step classical fourth order RK integrator
          if(ncase.eq.3)then
            call rkc4a(neqn,t0,tf,u0,nsteps,u)
          end if
C
C         Variable step classical fourth order RK integrator
          if(ncase.eq.4)then
            call rkc4b(neqn,t0,tf,u0,nsteps,abserr,relerr,u)
          end if
C
C         Fixed step RK Fehlberg (RKF45) integrator
          if(ncase.eq.5)then
            call rkf45a(neqn,t0,tf,u0,nsteps,u)
          end if
C
C         Variable step Fehlberg (RKF45) integrator
          if(ncase.eq.6)then
            call rkf45b(neqn,t0,tf,u0,nsteps,abserr,relerr,u)
          end if
C
C Advance solution
        t0=tf
        tf=tf+tp
        do i=1,neqn
          u0(i)=u(i)
        end do
C
C Next output
        end do
C
C Next integrator
      end do
C
C End of pdenonm
      end
