Name : MIRKDC

Authors: Wayne Enright (enright@cs.toronto.edu), 
         Paul Muir (muir@stmarys.ca).
Additional programming assistance: Mark Adams, Nicole DeYoung.
Thanks to Beta-Testers: Ian Gladwell, Larry Shampine, Richard Pancer.
 
Note: This code makes use of the Level 1 BLAS routines. These must
be linked with the MIRKDC package in order for it to work.


Purpose : Solves a system of first order, non-linear, boundary value, 
          ordinary differential equations, y'(t) = f(t,y(t)), with 
          separated boundary conditions, g_a(y(a))=0 and g_b(y(b))=0.
          Based on a mesh of points which subdivides [a,b], the code 
          uses mono-implicit Runge-Kutta (MIRK) methods to discretize 
          the ODE's and solves the resultant non-linear algebraic equations, 
          using a hybrid damped Newton/fixed Jacobian iteration, to get a
          discrete solution approximation at the meshpoints. Continuous 
          MIRK schemes are used to augment the discrete solution with a 
          C^1 continuous interpolant, u(t), which is then used to compute 
          an estimate of the solution defect, | u'(t) - f(t,u(t)) |. The 
          user defined tolerance is then applied to the defect as the 
          termination criterion.  Adaptive mesh redistribution is performed 
          to equidistribute the estimate of the defect. For further details 
          see [Enright,Muir, "Runge-Kutta Software with Defect Control for
          Boundary Value ODES",SIAM J. Sci. Comput. 17 (1996),479-497].


Usage   : call mirkdc(method, tol, neqns, Nsub, MxNsub, mesh_input, 
                      mesh, sol_input, Y, ldY, output_control, info, 
                      iwork, work, init_de, init_Y, fsub, gsub, dfsub, dgsub)


Arguments :
  
  method         - (input) - integer.
                   Possible values for `method' and the corresponding MIRK 
                   schemes are:
 
                   Value of `method' |   MIRK formula
                   ----------------------------------------------------------
                        2            |   Second order MIRK scheme
                                     |          (Trapezoidal scheme).
                        4            |   Fourth order MIRK scheme
                                     |          (Lobatto collocation scheme).
                        6            |   Sixth order MIRK scheme.
 
  tol            - (input) - double precision.
                   A user-defined tolerance applied to an estimate of the
                   defect of the approximate solution; the defect is the amount 
                   by which the continuous approximate solution fails to 
                   satisfy the ODE system. MIRKDC attempts to return, 
                   a continuous solution approximation, u(t), such that
                             |u'(t)-f(t,u(t))|/(|f(t,u(t))|+1) 
                   is less than `tol'. (This is a blended absolute/relative 
                   measure of the defect.) The same tolerance is applied to 
                   each component of the defect. 

  neqns          - (input) - integer.
                   The number of first order differential equations.
                   This is also the total number of boundary conditions.

  Nsub           - (input/output) - integer.
                   When passed to MIRKDC it is the initial user-defined number 
                   of subintervals into which the problem interval is 
                   partitioned.
                   Upon return from MIRKDC it is the number of subintervals 
                   in the final mesh. 

  MxNsub         - (input) - integer.
                   The user-defined maximum allowable number of subintervals. 
                   The value of this parameter will influence storage
                   requirements - see e.g. `iwork' and `work'.

  mesh_input     - (input) - integer.
                   a) A value of 0 indicates that a uniform initial 
                   mesh is to be set up by MIRKDC. 
                   b) A value of 1 indicates that an initial mesh is already 
                   available in the array 'mesh'. This option should be 
                   chosen if the user knows which regions of the interval 
                   correspond to difficult solution behavior, e.g. a boundary 
                   layer; the mesh should be finer in those regions.
                   c) A value of 2 indicates that the user is employing a 
                   ``continuation sequence'' of problems and that the mesh from
                   a previous problem will be available in the array `mesh'.   
                  
  mesh           - (input/output) - double precision array `(0:MxNsub)'.
                   On input to MIRKDC this is the set of points that partition 
                   the problem interval ( empty if `mesh_input' is 0 ). 
                   On successful return from MIRKDC it is the set of points 
                   defining the final mesh. 

  sol_input      - (input) - integer.
                   a) A value of 1 indicates that the initial solution estimate
                   will be provided by the user through the `init_Y' routine. 
                   This is the usual value for `sol_input' and the usual mode of
                   execution for MIRKDC.
                   It is recommended that `Y' not be set to a constant in case 
                   there is a singularity in f(t,y(t)), for that constant value.
                   b) A value of 2 indicates that the initial solution 
                   estimate at each mesh point is available in the 
                   array `Y'. This option should only be chosen when solving a 
                   continuation sequence of problems, where `Y' takes its value
                   from the solution of a previous problem.

  Y              - (input/output) - double precision two dimensional array of
                                                   size `ldY x (MxNsub+1)'.
                   On input to MIRKDC, if `sol_input'= 2, this is the discrete 
                   approximation to the solution, evaluated at each point of 
                   the mesh contained in `mesh'. (Empty if `sol_input'=1.)
                   The i-th column of `Y' contains the solution approximation 
                   vector for the i-th meshpoint.
                   On successful return from MIRKDC, this is the discrete
                   approximation to the solution evaluated at the points of the
                   final mesh. 

  ldY            - (input) - integer.
                   The leading dimension of the array `Y'. (`ldY' >= `neqns').

  output_control - (input) - integer.
                   Controls the amount of information output. Newton iteration 
                   counts, mesh selection information, and defect estimates are
                   examples of such information.
                   Possible values of `output_control' are :
                   0 - No output,  1 - Intermediate output,  2 - Full output.

  info           - (output) - integer. 
                   A communication flag:
                       `info' = 0 implies successful termination; the solution
                        has been obtained such that the estimates of the defect
                        satisfy the user defined tolerance - see `tol'.
                       `info' = 1 implies unsuccessful termination; the size 
                       of the new mesh would be too large. Unless `MxNsub' is
                       actually fairly small, this is usually an indication that
                       the code encountered difficulty with the convergence of
                       the Newton iterations. The user-supplied routines
                       `fsub', `dfsub', `gsub', and `dgsub' should be checked.

  iwork          - (output) - integer (work) array of size at least 
                                                         `neqns * (MxNsub + 1)'.

  work           - (output) - double precision (work) array of size at least  
                      `MxNsub * neqns * (2*neqns+35) + 6 * neqns * (4*neqns+1)'.

  init_de        - User-supplied subroutine to initialize the problem
                   dependent variables, 'leftbc', 'a', and 'b'.

                   call init_de(leftbc, a, b), where
                       leftbc - (output) - integer.
                                Number of boundary conditions applied at the
                                left endpoint of the problem interval.
                       a      - (output) - double precision.
                                The left endpoint of the problem interval.
                       b      - (output) - double precision.
                                The right endpoint of the problem interval.

  init_Y         - User-supplied subroutine to initialize the discrete 
                   solution approximation, `Y'. (Needed if `sol_input'=1).

                   call init_Y(Nsub, neqns, mesh, Y, ldY), where
                       Nsub   - (input) - integer.
                                The number of subintervals in the initial mesh.
                       neqns  - (input) - integer.
                                The number of differential equations in the
                                first order system.
                       mesh   - (input) - double precision array (0:Nsub).
                                Points making up the partition of the problem
                                interval.
                       Y      - (output) - double precision two dimensional
                                array of size `ldY x (Nsub + 1)'. 
                                The initial approximation to the discrete
                                solution at the points in the initial mesh.
                                The i-th column of `Y' is the solution 
                                approximation vector at mesh point, `mesh(i-1)'.
                       ldY    - (input) - integer.
                                The leading dimension of Y. (`ldY' >= `neqns').
                
  fsub           - User-supplied subroutine which defines f(t,y) for the 
                   first order system of differential equations, y' = f(t,y). 

                   call fsub(neqns, t, y, f), where
                       neqns   - (input) - integer.
                                 The number of differential equations in the
                                 first order system.
                       t       - (input) - double precision.
                                 A point in the problem interval. 
                       y       - (input) - double precision array `(1:neqns)'.
                                 The current solution approximation at t.
                       f       - (output) - double precision array `(1:neqns)'.
                                 The value of f(t,y).
 
  gsub           - User-supplied subroutine which defines the boundary
                   condition equations.

                   call gsub(neqns, ya, yb, bc), where
                       neqns    - (input) - integer.
                                  The number of differential equations in the
                                  first order system. Also the total number
                                  of boundary conditions.
                       ya       - (input) - double precision array `(1:neqns)'.
                                  Current solution approximation at the left 
                                  endpoint.
                       yb       - (input) - double precision array `(1:neqns)'.
                                  Current solution approximation at the right 
                                  endpoint.
                       bc       - (output) - double precision array `(1:neqns)'.
                                  Boundary condition equations evaluated at
                                  `ya' and `yb'. The first `leftbc' components
                                  of `bc' are `g_a(ya)'; the remaining
                                  `neqns-leftbc' components of `bc' are 
                                  `g_b(yb)'.
 
  dfsub          - User-supplied subroutine which defines the Jacobian, df/dy,
                   of the system of differential equations.  Since the array 
                   `JJ' has already been set to zero, it is only necessary to 
                   set the non-zero elements.

                   call dfsub(neqns, t, y, JJ), where
                       neqns    - (input) - integer.
                                  The number of differential equations in the
                                  first order system.
                       t        - (input) - double precision.
                                  A point in the problem interval. 
                       y        - (input) - double precision array `(1:neqns)'.
                                  The current solution approximation at t.
                       JJ       - (output) - double precision two dimensional 
                                  array of size `neqns x neqns'.
                                  Contains the Jacobian, df/dy.

  dgsub          - User-supplied subroutine which defines the Jacobian of
                   the boundary conditions, that is d(bc)/dy. The array `dbc' 
                   has been set to zero prior to the call to 'dgsub', so it
                   is only necessary for the user to set the non-zero elements.

                   call dgsub(neqns, ya, yb, dbc), where
                       neqns    - (input) - integer.
                                  The number of differential equations in the
                                  first order system. Also the total number
                                  of boundary conditions.
                       ya       - (input) - double precision array `(1:neqns)'.
                                  Current solution approximation at the left 
                                  endpoint.
                       yb       - (input) - double precision array `(1:neqns)'.
                                  Current solution approximation at the right 
                                  endpoint.
                       dbc      - (output) - double precision two dimensional 
                                  array of size `neqns x neqns'.
                                  Contains the Jacobian of the boundary 
                                  conditions.
 

Remarks:                  
1.   On output, both the `iwork' and `work' arrays contain information 
     that is required in order for the user to access the continuous 
     solution approximation using the  `sol_eval' routine.

2.   The code makes use of two external packages of routines.
     The COLROW package [Diaz,Fairweather,Keast, "COLROW and 
     ARECO:Fortran packages for solving certain almost block 
     diagonal linear systems by modified row and column 
     elimination",ACM Trans. Math. Soft. 9 (1983),376-380]
     is included with the MIRKDC source code. The double 
     precision version of the level 1 BLAS routines,
     [Lawson,Hanson,Kincaid,Krogh,"Basic linear algebra sub-programs 
     for Fortran usage",ACM Trans. Math. Soft. 5 (1979), 308-323],
     is available from the netlib software repository and must 
     be linked with the MIRKDC package.

3.   The usual technique for controlling the accuracy or quality of an
     appropriate solution is to estimate and control the global error.
     MIRKDC, in contrast, estimates the defect, u'(t) - f(t,u(t)) (where 
     u(t) is the continuous approximate solution) and applies the user 
     defined tolerance to control the defect.  See `tol' parameter 
     documentation and [Enright,Muir, "Runge-Kutta Software with Defect 
     Control for Boundary Value ODES",SIAM J. Sci. Comput. 17 (1996),479-497].
        
4.   When using continuation, it should be noted that the continuation
     sequence must be chosen so that the Newton iterations in MIRKDC
     do not fail since the absence of an appropriate `init_Y' routine
     will be problematic. (See e.g. [Ascher, Christiansen, Russell, 
     "Collocation software for boundary value ODEs", Math. Comp., 33 
     (1979), 659-679], for an example of how a BVODE code is used in 
     a parameter continuation sequence.)


     For further details see also the internal documentation for MIRKDC. 

-------------------------------------------------------------------------------

Name : Sol_eval

Purpose : This routine provides the solution approximation and its derivative 
          at a given point within the problem interval.


Usage : call sol_eval(t,z,z_prime,iwork,work)


Arguments :

  t             - (input) - double precision.
                  The evaluation point within the problem interval.

  z             - (output) - double precision array `(1:neqns)'.
                  The value of the interpolant at `t'.

  z_prime       - (output) - double precision array `(1:neqns)'.
                  The value of the derivative of the interpolant at `t'.

  iwork         - (input) - integer array of size at least 
                                                       `neqns * (MxNsub + 1)'.
                  This array contains integer information employed or set in
                  MIRKDC and needed by `sol_eval'.

   work         - (input) double precision array of size at least  
                    `MxNsub * neqns * (2*neqns+35) + 6 * neqns * (4*neqns+1)'.
                  This array contains double precision information employed or
                  set in MIRKDC and needed by `sol_eval'.


------------------------------------------------------------------------------
Sample Driver program:


c==============================================================================
c ***Model of main program to use MIRKDC, for Swirling Flow III test problem***
c       "Swirling Flow III", [Ascher, Mattheij, and Russell, 
c       "Numerical Solution of Boundary Value Problems for Ordinary 
c       Differential Equations", Classics in Applied Mathematics Series, 
c       Society for Industrial and Applied Mathematics, Philadelphia, 1995].
c==============================================================================
c
c***declaration of constants***
c
              integer neqns, MxNsub, ldY
              parameter(neqns=6,MxNsub=3000,ldY=10)
c
c  neqns       - The number of differential equations.
c  MxNsub      - The maximum number of subintervals.
c  ldY         - The leading dimension of the matrix Y.
c
c***declaration of remaining parameters to mirkdc***
c
             integer                method
             double precision       tol
             integer                Nsub
             integer                mesh_input
             double precision       mesh(0:MxNsub)
             integer                sol_input
             double precision       Y(ldY,(MxNsub+1))
             integer                output_control
             integer                info
             double precision       work(MxNsub*neqns*(2*neqns+35)
     +                                      +6*neqns*(4*neqns+1))
c            Make sure 'work' is at least
c            MxNsub*neqns*(2*neqns+35) + 6*neqns*(4*neqns+1)
c
             integer                iwork(neqns*(MxNsub+1))
c            Make sure 'iwork' is at least neqns*(MxNsub+1)
c
             double precision epsilon
             common /ode/ epsilon
c            Problem dependent parameter
c
c                  Value of 'method' |   MIRK formula
c                  -------------------------------------
c                       2            |   Second order MIRK scheme
c                       4            |   Fourth order MIRK scheme
c                       6            |   Sixth order MIRK scheme
c
c  tol           - A user-defined tolerance applied to the defect of the
c                  computed solution; the defect is the amount by which 
c                  the continuous approximate solution fails to satisfy 
c                  the ODE system. If 'mirkdc' returns successfully, 
c                  |def(t)|/(|f(t,y(t))|+1) will be less than 'tol', where 
c                  y'(t)=f(t,y(t)) is the ODE and def(t)= y'(t)-f(t,y(t)).
c                  The defect has one component per solution component.
c                  The same tolerance is applied to each component of the 
c                  defect.
c
c  Nsub          - On input to 'mirkdc', the number of subintervals into which 
c                  the problem interval is partitioned.
c                  On successful return from 'mirkdc', the number of sub - 
c                  intervals in the final mesh.
c
c  mesh_input    - a) A value of 0 indicates that a uniform initial 
c                     mesh is to be set up by 'mirkdc'.
c
c                  b) A value of 1 indicates that an initial mesh is 
c                     already available in the array 'mesh'. This option
c                     should be chosen if the user knows which regions of
c                     the interval are problematic so a finer mesh may be
c                     used on these trouble areas.
c
c                  c) A value of 2 indicates that the user is doing a 
c                     continuation problem in which case the mesh from
c                     a previous run is stored in the array 'mesh'.
c 
c  mesh          - On input to 'mirkdc', the set of points that partition the 
c                  problem interval, initially empty if 'mesh_input' = 0
c                  On successful return from 'mirkdc', the set of points
c                  defining the final mesh.
c
c  sol_input     - a) A value of 1 indicates that an initial solution
c                     estimate is to be provided by the user through
c                     the 'init_Y' routine. It is recommended that Y not 
c                     be set to a constant as there could be a singularity
c                     at that point in the differential equation. 
c
c                  b) A value of 2 indicates that an initial solution
c                     estimate at each mesh point is available in the 
c                     matrix 'Y'. This option should only be chosen when
c                     doing continuation problems where Y takes its
c                     value from the solution of a previous run.
c                         
c  Y              - On input to 'mirkdc' if 'sol_input' = 2, Y is a matrix 
c                   containing the discrete approximation to the solution, 
c                   evaluated at each point of the mesh contained in 'mesh'.
c                   On successful return from 'mirkdc', Y is a matrix 
c                   containing the discrete approximation to the solution 
c                   evaluated at the points in 'mesh'.
c                 
c
c  output_control - Controls the amount of information output. 
c                   Profiling of Newton iteration counts, mesh selection, 
c                   relative defect estimate are all examples of such 
c                   information.
c                   Possible values of output_control are:
c                            0 - No output.
c                            1 - Intermediate output.
c                            2 - Full output.
c
c  info           - Communication flag: 
c
c                   info = 0  successful termination - solution obtained
c                             to within user specified tolerance.
c
c                   info = -1 unsuccessful termination - the size
c                             of the new mesh would be too large.
c
             external init_de,init_Y,fsub,gsub,dfsub,dgsub 
c
c
c  See descriptions in the sample routines for further details.
c
c  init_de      - User - supplied subroutine to initialize the problem
c                 dependent variables, 'leftbc', 'a', and 'b'.
c                 As well this is a good place to set up any common blocks 
c                 needed to pass ODE parameter related info to any other 
c                 user defined routines.
c                 call init_de(leftbc, a, b) where
c                 leftbc - Number of boundary conditions supplied at the
c                          left endpoint of the problem interval.(output)
c                 a      - The left endpoint of the interval.(output)
c                 b      - The right endpoint of the interval.(output)  
c
c  init_Y       - User - supplied subroutine to initialize the discrete 
c                 solution approximation, Y. Y is a matrix of size 
c                 (neqns) x (Nsub + 1). The i'th column of Y is the 
c                 solution approximation vector for the mesh point,
c                 mesh(i-1).
c                 call init_Y(Nsub, neqns, mesh, Y, ldY) where
c                 Nsub   - The number of subintervals in the mesh.(input)
c                 neqns  - The number of differential equations in the
c                          first order system.(input)
c                 mesh   - Points making up the partition of the problem
c                          interval.(input)
c                 Y      - The initial approximation of the discrete
c                          solution at the points in the initial mesh.
c                          (output).
c                 ldY    - The leading dimension of Y.(input)
c               
c  fsub         - User - supplied subroutine which defines f(t,y) for the 
c                 first order system of differential equations, y' = f(t,y). 
c                 call fsub(neqns, t, y, f) where
c                 neqns   - The number of differential equations in the
c                           first order system.(input)
c                 t       - A point in the problem interval.(input) 
c                 y       - The solution approximation at time, t.(input)
c                 f       - The value of f(t,y).(output)
c
c  gsub         - User - supplied subroutine which defines the boundary
c                 condition equations.
c                 call gsub(neqns, ya, yb, bc) where
c                 neqns    - The number of differential equations in the
c                            first order system.(input)
c                 ya       - Denotes the boundary conditions that have 
c                            been specified at the left endpoint.(input)
c                 yb       - Denotes the boundary condidtions that have
c                            been specified at the right endpoint.(input)
c                 bc       - The value of the boundary condition equation.
c                            (output)
c
c  dfsub        - User - supplied subroutine which defines the Jacobian of
c                 the system of differential equations with respect to y,
c                 i.e. df/dy. The Jacobian has already been set to zero so
c                 it is only necessary that the user input the non-zero
c                 elements.
c                 call dfsub(neqns, t, y, JJ) where
c                 neqns    - The number of differential equations in the 
c                            first order system.(input)
c                 t        - A point in the problem interval.(input)
c                 y        - The solution approximation at time, t.(input)
c                 JJ       - A two-dimensional array containing the
c                            Jacobian.(output)
c
c  dgsub        - User - supplied subroutine which defines the Jacobian of
c                 the boundary conditions, that is d bc/dy. The Jacobian 
c                 has been set to zero prior to the call to 'dgsub', so it
c                 is only necessary that the user input the non-zero
c                 elements.
c                 call dgsub(neqns, ya, yb, dbc) where
c                 neqns    - The number of differential equations in the 
c                            first order system.(input)
c                 ya       - Denotes the boundary conditions that have been
c                            specified at the left endpoint.(input)
c                 yb       - Denotes the boundary conditions that have been
c                            specified at the right endpoint.(input)
c                 dbc      - A two-dimensional array containing the Jacobian
c                            of the boundary conditions.(output)
c
c     ***declaration of local variables***
c
             integer                i,j,fail
             double precision       x,h
             double precision       z(neqns),zp(neqns)
c
c  i,j       - loop indexes
c
c  fail      - flag to control execution of 'sol_eval'
c
c   x        - local variable to store a meshpoint
c
c   h        - local variable to store the size of a subinterval 
c
c   z        - value of the interpolant evaluated at a mesh point
c
c   zp       - value of the interpolant's derivative evaluated at a 
c              mesh point
c-------------------------------------------------------------------------
c***Setup parameters***
c
      write(*,*)
      write(*,*)'Initialization :'
      write(*,*)  
c
c     Set the order of the MIRK scheme.
      write(*,*)
      write(*,*)'Input value of method 2, 4 or 6' 
      read(5,*) method
c
c     Set the defect tolerance. 
      write(*,*)'Input value of tolerance'
      read(5,*)tol
c
c     Set the number of subintervals in the initial mesh.
      write(6,*)'Input the initial number of subintervals.'
      Read(5,*) Nsub
c
c     Set the value of problem dependent parameter, epsilon
      write(6,*)'Input the value of epsilon'
      read(5,*)epsilon
c
c     Set mesh_input to indicate that mirkdc should begin with a uniform mesh.
      mesh_input = 0
c
c     Set indicator for type of initial solution estimate.
      sol_input = 1
c
c     Set value for output_control
      output_control = 1
c
      write(6,*)
      write(6,*)'Swirling Flow III'
      write(6,*)
c
c     ******Call mirkdc******
c      
      call mirkdc(method,tol,neqns,Nsub,MxNsub,mesh_input,
     +        mesh,sol_input,Y,ldY,output_control,info,iwork,work,
     +        init_de,init_Y,fsub,gsub,dfsub,dgsub)
c      
c     ******Return from mirkdc******
c
c     Upon successful return from mirkdc, print out solution.
c     Solution interpolant is available through the sol_eval routine        
c     This routine requires access to the work array.
c
c     Check return from MIRKDC
c
      if (info .NE. 0) then
                write(6,*)'Unsuccessful termination'
                stop
      endif
c
c     Compute solution approximations at 11 points.
c      
      h = (mesh(Nsub) - mesh(0))/10.0d0
      do 40 i = 1, 11
         x = (mesh(0) + (i-1)*h)
         call sol_eval(x,z,zp,fail,iwork,work)
         if (fail.EQ.0) then
          write(6,*)'At meshpoint x=', x,' the solution is'
          do 35 j=1,neqns
             write(6,100)z(j)
 100         format(1x,d30.15) 
 35       continue
         else
          write(6,*)'Attempted to evaluate the interpolant at a point', 
     +      'that was outside of the problem interval. No information',
     +      'is available for that point.'
         end if
 40   continue
 
      stop
      end
 
c A Test Problem for Mirkdc
c==============================================================================
c CONTENTS
c
c       This file contains the routines associated  with the nonlinear 
c       problem "Swirling Flow III", [Ascher, Mattheij, and Russell, 
c       "Numerical Solution of Boundary Value Problems for Ordinary 
c       Differential Equations", Classics in Applied Mathematics Series, 
c       Society for Industrial and Applied Mathematics, Philadelphia, 1995].
c       This is the special case of counter rotating disks - thus Sigma_0=-1
c       and Sigma_1=1. It is stated as:       
c                         
c               Original Form                First Order System Form
c               -------------                -----------------------
c                                             y1' = y2,
c               g'' = gf' - fg'                                
c                     --------- ,             y2' = y1*y4 - y3*y2
c                      epsilon                      ------------- ,
c                                                      epsilon
c               f'''' = -ff''' - gg'          y3' = y4 ,
c                       ------------ ,
c                         epsilon             y4' = y5 ,  y5' = y6 ,
c
c                                             y6' = -y3*y6 - y1*y2
c                                                   -------------- ,
c                                                       epsilon
c       with
c
c               g(0) = -1, g(1) = 1,             y1(0) = -1, y1(1) = 1,
c
c       f(0) = f'(0) = f(1) = f'(1) = 0,         y3(0) = y4(0) = 0,
c                                                y3(1) = y4(1) = 0,
c       and             
c             epsilon = {a small parameter}.
c
c       No closed form solution is available.
c==============================================================================
c
        subroutine init_de(leftbc,a,b)

c	The purpose of this routine is to initialize the problem dependent
c	variables 'leftbc', the number of boundary conditions at the left
c	endpoint of the problem interval, 'neqns', the number of differential 
c	equations, 'a', the left endpoint of the problem interval, and 'b',
c	the right endpoint of the problem interval.
c------------------------------------------------------------------------------
c
c***declaration of parameters***
c	    exports:
		integer leftbc
		double precision a,b

c		'leftbc' the number of boundary conditions imposed at the
c			 left end of the problem interval
c		'a' the left endpoint of the problem interval
c		'b' the right endpoint of the problem intervalc
c------------------------------------------------------------------------------
c
        leftbc = 3
c
        a = 0.0d0
        b = 1.0d0
c
        return
        end
c       
c=============================================================================
c
        subroutine init_Y(Nsub,neqns,mesh,Y,ldY)
c
c	The purpose of this routine is to initialize the discrete solution 
c	approximation, 'Y'. The input, 'Nsub', is the number of subintervals 
c	in the initial partition of the problem interval. The points of the
c	partition are stored in 'mesh'. 'Y' has 'neqns' components per meshpoint
c-----------------------------------------------------------------------------
c		
c***declaration of parameters***
c	    imports:
		integer Nsub, neqns, ldY
		double precision mesh(0:Nsub)
c
c		'Nsub' number of subintervals in 'mesh'
c		'neqns' number of differential equations
c		'mesh' points making up partition of problem interval
c               'ldY' the leading dimension of Y.
c
c	    exports:
		double precision Y(ldY,(Nsub+1))
c
c		'Y' the initial approximation to the discrete solution
c		    at the points of the initial mesh.
c
c***declaration of local variables***
c
		integer i
c
c		'i' loop index from 0 to Nsub
c------------------------------------------------------------------------------
c       The initial approximation is based on straight lines joining the
c       boundary conditions for each differential equation in the first 
c       order system. Boundary conditions are given for the first,
c       third and fourth equations, defining the lines y = 2t-1 (or mesh(i))
c       y = 0 and y = 0. Since there are no boundary conditions for the second,
c       fifth, or sixth equations, the line drawn for y2 is y = 1 (the
c       derivative of y = t since y1' = y2 ) and y5 and y6 are initialized
c       by the line y = 0.
c----------------------------------------------------------------------------
c  
        do 10 i = 0, Nsub
            Y(1, i+1) = 2.0d0*mesh(i)-1.0d0
            Y(2, i+1) = 2.d0
            Y(3, i+1) = 0.d0
            Y(4, i+1) = 0.d0
            Y(5, i+1) = 0.d0
            Y(6, i+1) = 0.d0   
 10     continue
c
        return
        end
c
c=============================================================================
c
        subroutine fsub(neqns,t,y,f)
c
c       This routine defines f(t,y) for the first order system of 
c       differential equations.`neqns' is the size of the system.
c-----------------------------------------------------------------------------
c
c***declaration of parameters***
c            imports:
                integer                    neqns
                double precision           t, y(neqns)
c
c               `neqns' Number of differential equations.
c               `t' A point in the problem interval where f(t,y(t)) is to be
c                   evaluated.
c               `y' Current solution approximation at `t', used in the 
c                   evaluation of f(t,y(t)).
c            exports:
                double precision           f(neqns)
c
c               `f' Value of f(t,y(t)) for given `t' and `y'.
c
c***declaration for parameter 'epsilon'***
                double precision epsilon
                common /ode/ epsilon
c------------------------------------------------------------------------------
c        
        f(1) = y(2)
        f(2) = (y(1)*y(4) - y(3)*y(2))/epsilon 
        f(3) = y(4)
        f(4) = y(5)
        f(5) = y(6)
        f(6) = (-y(3)*y(6) -  y(1)*y(2))/epsilon
c
        return
        end
c
c==============================================================================
c
        subroutine dfsub(neqns,t,y,JJ)
c
c       This routine defines the Jacobian, d f(t,y(t))/dy, of the system of 
c       differential equations. Since `JJ' has already been set to zero, it is 
c       only necessary to set the non-zero elements.
c------------------------------------------------------------------------------
c
c***declaration of parameters***
c            imports:
                integer                neqns
                double precision       t, y(neqns)
c
c               `neqns' Number of differential equations.
c               `t' A point in the problem interval where d f(t,y(t))/dy is to 
c                   be evaluated.
c               `y' Current solution approximation at `t', used in the 
c                   evaluation of d f(t,y(t))/dy.
c            exports:
                double precision       JJ(neqns,neqns)
c
c               `JJ' Value of d f(t,y(t))/dy for given `t' and `y'.
c
c       declaration for parameter 'epsilon'
                double precision epsilon
                common /ode/ epsilon
c------------------------------------------------------------------------------
c
        JJ(1,2) = 1.0d0
c
        JJ(2,1) = y(4)/epsilon
        JJ(2,2) = -y(3)/epsilon
        JJ(2,3) = -y(2)/epsilon
        JJ(2,4) = y(1)/epsilon
c       
        JJ(3,4) = 1.0d0
c       
        JJ(4,5) = 1.0d0
c        
        JJ(5,6) = 1.0d0
c
        JJ(6,1) = -y(2)/epsilon
        JJ(6,2) = -y(1)/epsilon
        JJ(6,3) = -y(6)/epsilon
        JJ(6,6) = -y(3)/epsilon
c
        return
        end
c
c==============================================================================
c
        subroutine gsub(neqns,ya,yb,bc)
c
c       This routine  computes  the 'neqns' boundary conditions, assumed
c	to be separated, applied at either the left endpoint, a, or the
c	right endpoint, b.  Each boundary condition is assigned to the
c	corresponding component of 'bc' as follows. The code will attempt
c	to find 'ya' and 'yb' values (i.e. ya = y(a) and yb = y(b)) in order
c	to make each component of 'bc' equal to zero. This each boundary 
c	condition must be written as some quantity set equal to 0. Thus,
c	the boundary condition, ya(2) = 3, would be assigned to say the
c	first boundary condition as, bc(1) = ya(2) - 3.
c------------------------------------------------------------------------------
c
c***declaration of parameters***
c            imports:
                integer                   neqns
                double precision          ya(neqns),yb(neqns)
c
c               `neqns' Number of differential equations.
c               `ya' Current solution approximation at `a', used in the 
c                   evaluation of g_a(y(a)).
c               `yb' Current solution approximation at `b', used in the 
c                   evaluation of g_b(y(b)).
c            exports:
                 double precision          bc(neqns)
c
c                `bc' The first `leftbc' locations contain g_a(y(a)).
c                     The remaining `neqns-leftbc' locations contain g_b(y(b)).
c------------------------------------------------------------------------------
c
        bc(1) = ya(1) + 1.0d0
        bc(2) = ya(3)
        bc(3) = ya(4)
c
        bc(4) = yb(1) - 1.0d0
        bc(5) = yb(3)
        bc(6) = yb(4)
c
        return
        end
c
c==============================================================================
c
        subroutine dgsub(neqns,ya,yb,dbc)
c
c       This routine defines the Jacobian, d bc/dy, of the boundary conditions.
c       `dbc' has already been set to zero so it is only necessary to set
c       the non-zero elements.
c------------------------------------------------------------------------------
c
c***declaration of parameters***
c            imports:
                integer                   neqns
                double precision          ya(neqns),yb(neqns)
c
c               `neqns' Number of differential equations.
c               `ya' Current solution approximation at `a', used in the 
c                   evaluation of d g_a(y(a))/d ya.
c               `yb' Current solution approximation at `b', used in the 
c                   evaluation of d g_b(y(b))/ d yb.
c            exports:
                double precision          dbc(neqns,neqns)
c
c               `dbc' Value of d bc/dy for given `t',`ya', and `yb'.
c---------------------------------------------------------------------------
c
        dbc(1,1) = 1.0d0
        dbc(2,3) = 1.0d0
        dbc(3,4) = 1.0d0
        dbc(4,1) = 1.0d0
        dbc(5,3) = 1.0d0
        dbc(6,4) = 1.0d0
c
        return
        end