      subroutine intpar(neqn,nout,nsteps,t0,tf,abserr,relerr)
C
C Subroutine intpar sets the parameters to control the
C integration of the 2 x 2 ODE system
C
C Double precision coding is used 
      implicit double precision(a-h,o-z)
C
C Number of ODEs
      neqn=2
C
C Number of output points
      nout=6
C
C Maximum number of steps in the interval t0 to tf
      nsteps=100
C
C Initial, final values of the independent variable
      t0=0.0d0
      tf=1.0d0
C
C Error tolerances
      abserr=1.0d-05
      relerr=1.0d-05
      return
C
C End of intpar
      end
      subroutine inital(neqn,t,u0)
C
C Subroutine inital sets the initial condition vector
C for the 2 x 2 ODE 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 Initial condition
      u0(1)=0.0d0
      u0(2)=2.0d0
      return
C
C End of inital
      end
      subroutine par(a,b)
C
C Subroutine par sets the parameters for the 2 x 2 ODE
C problem
C
C Double precision coding is used
      implicit double precision(a-h,o-z)
C
C Problem parameters
      a=5.5d0
      b=4.5d0
      return
C
C End of par
      end
      subroutine derv(neqn,t,u,ut)
C
C Subroutine derv computes the derivative vector
C of the 2 x 2 ODE problem
C
C Double precision coding is used
      implicit double precision(a-h,o-z)
C
C Size the arrays
      dimension u(neqn), ut(neqn)
C
C Problem parameters
      call par(a,b)
C
C Derivative vector
      ut(1)=-a*u(1)+b*u(2)
      ut(2)= b*u(1)-a*u(2)
      return
C
C End of derv
      end
      subroutine fprint(no,ncase,neqn,t,u)
C
C Subroutine fprint displays the numerical and
C analytical solutions to the 2 x 2 ODE problem
C
C Double precision coding is used
      implicit double precision(a-h,o-z)
C
C Size the arrays
      dimension u(neqn)
C
C Problem parameters
      call par(a,b)
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(/,6x,'euler2a integrator')
C
C         Variable step modified Euler
          else if(ncase.eq.2)then
            write(no,12)
 12         format(/,6x,'euler2b integrator')
C
C         Fixed step classical fourth order RK
          else if(ncase.eq.3)then
            write(no,13)
 13         format(/,6x,'rkc4a integrator')
C
C         Variable step classical fourth order RK
          else if(ncase.eq.4)then
            write(no,14)
 14         format(/,6x,'rkc4b integrator')
C
C         Fixed step RK Fehlberg 45
          else if(ncase.eq.5)then
            write(no,15)
 15         format(/,6x,'rkf45a integrator')
C
C         Variable step RK Fehlberg 45
          else if(ncase.eq.6)then
            write(no,16)
 16         format(/,6x,'rkf45b integrator')
          end if
C
C Heading
      write(no,1)ncase
 1    format(/,5x,' int = ',i3)
      write(no,2)
 2    format(/,9x,'t',3x,'u1(num)',4x,'u1(ex)',8x,'diff1',/,
     1        10x,    3x,'u2(num)',4x,'u2(ex)',8x,'diff2',/)
C
C End of t = 0 heading
      end if
C
C Analytical solution
      u1exact=dexp(-(a-b)*t)-dexp(-(a+b)*t)
      u2exact=dexp(-(a-b)*t)+dexp(-(a+b)*t)
C
C Difference between exact and numerical solution vectors
      diff1=u(1)-u1exact
      diff2=u(2)-u2exact       
C
C Display the numerical and exact solutions, and their difference
      write(no,3)t,u(1),u1exact,diff1,u(2),u2exact,diff2
 3    format(f10.2,2f10.5,e13.4,/,10x,2f10.5,e13.4,/)
      return
C
C End of fprint
      end
