      subroutine intpar(neqn,nout,nsteps,t0,tf,abserr,relerr)
C
C Subroutine intpar sets the parameters to control the
C integration of the nonlinear PDE problem
C
C Double precision coding is used 
      implicit double precision(a-h,o-z)
C
C Number of ODEs
      neqn=22
C
C Number of output points
      nout=6
C
C Maximum number of steps in the interval t0 to tf
      nsteps=1000
C
C Initial, final values of the independent variable
      t0=0.0d0
      tf=0.2d0
C
C Error tolerances
      abserr=1.0d-01
      relerr=1.0d-05
      return
C
C End of intpar
      end
      subroutine par
C
C Double precision coding is used
      implicit double precision(a-h,o-z)
C
C Global variables
      common/p/  alpha,     xk,  sigma,     xl,     ui,     ua,
     +               a,      e,   rhos,    cps,    xls,     cs
C
C Subroutine par sets the parameters for the WTC model
      alpha=1.0d-06
      xk=1.0d0
      sigma=5.67d-08  
      xl=0.1d0
      ui=298.0d0
      ua=2000.0d0
      a=1.0d0
      e=1.0d0
      rhos=7800.0d0
      cps=435.0d0
      xls=0.025d0
      cs=rhos*cps*xls
      return
C
C End of par
      end
      subroutine inital(neqn,t,u0)
C
C Subroutine inital sets the initial condition vector 
C for the nonlinear PDE problem
C
C Double precision coding is used
      implicit double precision(a-h,o-z)
C
C Global variables
      common/p/  alpha,     xk,  sigma,     xl,     ui,     ua,
     +               a,      e,   rhos,    cps,    xls,     cs
C
C Size the arrays
      dimension u0(neqn)
C
C Problem parameters
      call par
C
C Initial condition
      do i=1,neqn
        u0(i)=ui
      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 nonlinear PDE problem
C
C Double precision coding is used
      implicit double precision(a-h,o-z)
C
C Global variables
      common/p/  alpha,     xk,  sigma,     xl,     ui,     ua,
     +               a,      e,   rhos,    cps,    xls,     cs
C
C Size the arrays
      dimension u(neqn), ut(neqn)
C
C Spatial grid
      dx=1.0d0/dfloat(neqn-2)
      dxs=dx*dx
C
C Insulation
      do i=1,neqn-2
        if(i.eq.1)then
          u0=u(2)+2.0d0*dx*xl*(sigma/xk)*(a*ua**4-e*u(1)**4)
          ut(i)=(u(2)-2.0d0*u(1)+u0)/dxs
        else 
          ut(i)=(u(i+1)-2.0d0*u(i)+u(i-1))/dxs
        end if
      end do
C
C Steel
      ut(neqn)=(xk*xl/alpha)*(1.0d0/cs)*(u(neqn-2)-u(neqn-1))/dx
      ut(neqn-1)=ut(neqn)
C
C End of derv
      return
      end
      subroutine fprint(no,ncase,neqn,t,u)
C
C Subroutine fprint displays the numerical solution to the 
C nonlinear PDE problem
C
C Double precision coding is used
      implicit double precision(a-h,o-z)
C
C Global variables
      common/p/  alpha,     xk,  sigma,     xl,     ui,     ua,
     +               a,      e,   rhos,    cps,    xls,     cs
C
C Size the arrays
      dimension u(neqn)
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,2)ncase,neqn
 2    format(/,' ncase = ',i2,'   neqn = ',i2,/)
      write(no,3)
 3    format('     t      u(1)     u(im)    u(neqn-1)')
C
C End of t = 0 heading
      end if
C
C Numerical and analytical solution output 
C  	
C     Grid index of midpoint 
      im=neqn/2
C	
C     Display the numerical solution
      tmin=t*xl**2/alpha/60.0d0
      write(no,4)  tmin,     u(1),    u(im),   u(neqn-1)
 4    format(f6.1,3f10.2)
      return
C
C End of fprint
      end
