c
c----------------------------------------------------------------
c
c  problem 1 - see companion paper
c
      real fspace(2000), zeta(4), tol(2), z(4), u(4), err(4)
      integer ispace(200), m(1), ipar(11), ltol(2)
      external fsub, dfsub, gsub, dgsub
c
      write (6,99999)
c
c     one differential equation of order 4.
      m(1) = 4
c     give location of boundary conditions
      zeta(1) = 1.
      zeta(2) = 1.
      zeta(3) = 2.
      zeta(4) = 2.
c     set up parameter array.
c     use default values for all parameters except for initial
c     mesh size, no. of tolerances and sizes of work arrays
      do 10 i=1,11
        ipar(i) = 0
   10 continue
      ipar(3) = 1
      ipar(4) = 2
      ipar(5) = 2000
      ipar(6) = 200
c     two error tolerances (on u and its second derivative)
      ltol(1) = 1
      ltol(2) = 3
      tol(1) = 1.e-7
      tol(2) = 1.e-7
c
      call colsys(1, m, 1., 2., zeta, ipar, ltol, tol, dummy, ispace,
     * fspace, iflag, fsub, dfsub, gsub, dgsub, dummy)
c
      if (iflag.ne.1) stop
c     calculate the error at 101 points using the known
c     exact solution
      x = 1.
      do 20 i=1,4
        err(i) = 0.
   20 continue
      do 40 j=1,101
        call appsln(x, z, fspace, ispace)
        call exact(x, u)
        do 30 i=1,4
          err(i) = amax1(err(i),abs(u(i)-z(i)))
   30   continue
        x = x + .01
   40 continue
      write (6,99998) (err(i),i=1,4)
      stop
99999 format (1h1, 35h example of a simple problem setup./10h  uniforml,
     * 36hy loaded beam of variable stiffness,/21h  simply supported at,
     * 11h both ends./)
99998 format (/27h error tolerances satisfied//22h the exact errors are,
     * /7x, 4e12.4)
      end
c................................................................
      subroutine fsub(x, z, f)
      real z(4), f(1)
      f(1) = (1.-6.*x**2*z(4)-6.*x*z(3))/x**3
      return
      end
c................................................................
      subroutine dfsub(x, z, df)
      real z(4), df(1,4)
      df(1,1) = 0.
      df(1,2) = 0.
      df(1,3) = -6./x**2
      df(1,4) = -6./x
      return
      end
c................................................................
      subroutine gsub(i, z, g)
      real z(4)
      go to (10, 20, 10, 20), i
   10 g = z(1) - 0.
      return
   20 g = z(3) - 0.
      return
      end
c................................................................
      subroutine dgsub(i, z, dg)
      real z(4), dg(4)
      do 10 j=1,4
        dg(j) = 0.
   10 continue
      go to (20, 30, 20, 30), i
   20 dg(1) = 1.
      return
   30 dg(3) = 1.
      return
      end
c................................................................
      subroutine exact(x, u)
      real u(4)
c     exact solution
      u(1) = .25*(10.*alog(2.)-3.)*(1.-x) + .5*(1./x+(3.+x)*alog(x)-x)
      u(2) = -.25*(10.*alog(2.)-3.) + .5*(-1./x/x+alog(x)+(3.+x)/x-1.)
      u(3) = .5*(2./x**3+1./x-3./x/x)
      u(4) = .5*(-6./x**4-1./x/x+6./x**3)
      return
      end
c
c----------------------------------------------------------------
c
c  problem 2 - see companion paper
c
      real zeta(4), fspace(40000), tol(4), z(4)
      integer m(2), ipar(11), ispace(2500), ltol(4)
      common eps, dmu, eps4mu, gamma, xt
      external solutn, fsub, dfsub, gsub, dgsub
c     define constants, print a heading.
      gamma = 1.1
      eps = .01
      dmu = eps
      eps4mu = eps**4/dmu
      xt = sqrt(2.*(gamma-1.)/gamma)
      write (6,99999) gamma, xt, eps, dmu, eps4mu
c     define no. of differential equations.
      ncomp = 2
c     orders
      m(1) = 2
      m(2) = 2
c     interval ends
      aleft = 0.
      aright = 1.
c     locations of side conditions
      zeta(1) = 0.
      zeta(2) = 0.
      zeta(3) = 1.
      zeta(4) = 1.
c     ipar  values
c     a nonlinear problem
      ipar(1) = 1
c     4 collocation points per subinterval
      ipar(2) = 4
c     initial uniform mesh of 10 subintervals
      ipar(3) = 10
      ipar(8) = 0
c     dimension of real work array  fspace  is 40000
      ipar(5) = 40000
c     dimension of integer work array  ispace  is 2500
      ipar(6) = 2500
c     (these dimensions of  fspace  and  ispace
c      enable  colsys  to use meshes of up to 192 intervals.)
c     print full output.
      ipar(7) = -1
c     initial approximation for nonlinear iteration is provided
c     in  solutn
      ipar(9) = 1
c     a regular problem
      ipar(10) = 0
c     no fixed points in the mesh
      ipar(11) = 0
c     tolerances on  all components
      ipar(4) = 4
      do 10 i=1,4
        ltol(i) = i
        tol(i) = 1.e-5
   10 continue
c     call  colsys
      call colsys(ncomp, m, aleft, aright, zeta, ipar, ltol, tol,
     * fixpnt, ispace, fspace, iflag, fsub, dfsub, gsub, dgsub, solutn)
c     print values of the obtained approximate solution at points
c     x = 0,.05, ..., 1.
      x = 0.
      write (6,99998)
      np1 = 21
      do 20 iii=1,np1
        call appsln(x, z, fspace, ispace)
        write (6,99997) x, z
        x = x + .05
   20 continue
      stop
99999 format (1h1, 27hdimpling of spherical caps./8h gamma =,
     * f7.2/6h  xt =, e12.5/6h eps =, e12.5/6h  mu =, e12.5/9h eps**4/m,
     * 3hu =, e12.5)
99998 format (1h1, 44h       x             phi           dphi     ,
     * 23h      psi          dpsi/)
99997 format (6x, f5.2, 4x, 6e15.5)
      end
c................................................................
      subroutine solutn(x, z, dmval)
      common eps, dmu, eps4mu, gamma, xt
      dimension z(4), dmval(2)
      cons = gamma*x*(1.-.5*x*x)
      dcons = gamma*(1.-1.5*x*x)
      d2cons = -3.*gamma*x
      if (x.gt.xt) go to 10
      z(1) = 2.*x
      z(2) = 2.
      z(3) = -2.*x + cons
      z(4) = -2. + dcons
      dmval(2) = d2cons
      go to 20
   10 z(1) = 0.
      z(2) = 0.
      z(3) = -cons
      z(4) = -dcons
      dmval(2) = -d2cons
   20 dmval(1) = 0.
      return
      end
c................................................................
      subroutine fsub(x, z, f)
      dimension z(4), f(2)
      common eps, dmu, eps4mu, gamma, xt
      f(1) = z(1)/x/x - z(2)/x + (z(1)-z(3)*(1.-z(1)/x)-gamma*x*(1.-x*x/
     * 2.))/eps4mu
      f(2) = z(3)/x/x - z(4)/x + z(1)*(1.-z(1)/2./x)/dmu
      return
      end
c................................................................
      subroutine dfsub(x, z, df)
      dimension z(4), df(2,4)
      common eps, dmu, eps4mu, gamma, xt
      df(1,1) = 1./x/x + (1.+z(3)/x)/eps4mu
      df(1,2) = -1./x
      df(1,3) = -(1.-z(1)/x)/eps4mu
      df(1,4) = 0.
      df(2,1) = (1.-z(1)/x)/dmu
      df(2,2) = 0.
      df(2,3) = 1./x/x
      df(2,4) = -1./x
      return
      end
c................................................................
      subroutine gsub(i, z, g)
      dimension z(4)
      go to (10, 20, 10, 30), i
   10 g = z(1)
      return
   20 g = z(3)
      return
   30 g = z(4) - .3*z(3) + .7
      return
      end
c................................................................
      subroutine dgsub(i, z, dg)
      dimension z(4), dg(4)
      do 10 j=1,4
        dg(j) = 0.
   10 continue
      go to (20, 30, 20, 40), i
   20 dg(1) = 1.
      return
   30 dg(3) = 1.
      return
   40 dg(4) = 1.
      dg(3) = -.3
      return
      end
c
c----------------------------------------------------------------
c
c  problem 3 - see companion paper
c
      real zeta(5), fspace(40000), tol(2), sval(3), elval(3)
      integer m(2), ispace(2500), ltol(2), ipar(11)
      real z(5)
      common en, s, el, cons
      external fsub, dfsub, gsub, dgsub, solutn
      data sval /.2,.1,.05/, elval /60.,120.,200./
c
      en = .2
      cons = .5*(3.-en)
      ncomp = 2
      m(1) = 2
      m(2) = 3
      aleft = 0.
      aright = 1.
c
      zeta(1) = 0.
      zeta(2) = 0.
      zeta(3) = 0.
      zeta(4) = 1.
      zeta(5) = 1.
c
      ipar(1) = 1
      ipar(2) = 4
      ipar(3) = 10
      ipar(4) = 2
      ipar(5) = 40000
      ipar(6) = 2500
      ipar(7) = 0
      ipar(8) = 0
      ipar(9) = 1
      ipar(10) = 0
      ipar(11) = 0
c
      ltol(1) = 1
      ltol(2) = 3
      tol(1) = 1.e-5
      tol(2) = 1.e-5
c
c     solve a chain of 3 problems
      do 30 ijk=1,3
        s = sval(ijk)
        el = elval(ijk)
        if (ijk.eq.1) go to 10
c        set continuation parameters
        ipar(9) = 3
        ipar(3) = ispace(1)
   10   continue
        write (6,99999) en, s, el
c
        call colsys(ncomp, m, aleft, aright, zeta, ipar, ltol, tol,
     *   fixpnt, ispace, fspace, iflag, fsub, dfsub, gsub, dgsub,
     *   solutn)
c
        if (iflag.ne.1) stop
c        print values of the obtained approximate solution at points
c        x = 0,1,2, ..., l.
        is6 = ispace(6)
        is5 = ispace(1) + 2
        x = 0.
        write (6,99998)
        np1 = el + 1.5
        do 20 iii=1,np1
          call approx(ii, x, z, fspace(is6), fspace(1), ispace(1),
     *     fspace(is5), ispace(2), ncomp, m, ispace(4), 1, dm, 0)
          xl = x*el
          z(2) = z(2)/el
          z(4) = z(4)/el
          z(5) = z(5)/el/el
          write (6,99997) xl, z
          x = x + 1./el
   20   continue
   30 continue
      stop
99999 format (1h1, 38h rotating flow over a stationary disk./8h  parame,
     * 11hters -  n =, f5.2, 6h   s =, f5.2, 6h   l =, f6.1/)
99998 format (1h1, 44h       x              g              dg     ,
     * 38h       h             dh            d2h/)
99997 format (6e15.5)
      end
c................................................................
      subroutine solutn(x, z, dmval)
      common en, s, el, cons
      real z(5), dmval(2)
      ex = exp(-el*x)
      z(1) = 1. - ex
      z(2) = el*ex
      z(3) = -el**2*x**2*ex
      z(4) = (el**3*x**2-2.*el**2*x)*ex
      z(5) = (-el**4*x**2+4.*el**3*x-2.*el**2)*ex
      dmval(1) = -el*z(2)
      dmval(2) = (el**5*x*x-6.*el**4*x+6.*el**3)*ex
      return
      end
c................................................................
      subroutine fsub(x, z, f)
      real z(1), f(1)
      common en, s, el, cons
      f(1) = -el*(cons*z(3)*z(2)+(en-1.)*z(4)*z(1)) + el**2*s*(z(1)-1.)
      f(2) = -el*(cons*z(3)*z(5)+en*z(4)**2) + el**2*s*z(4) +
     * el**3*(1.-z(1)**2)
      return
      end
c................................................................
      subroutine dfsub(x, z, df)
      real z(1), df(2,1)
      common en, s, el, cons
      df(1,1) = -el*(en-1.)*z(4) + el**2*s
      df(1,2) = -el*cons*z(3)
      df(1,3) = -el*cons*z(2)
      df(1,4) = -el*(en-1.)*z(1)
      df(1,5) = 0.
      df(2,1) = -el**3*2.*z(1)
      df(2,2) = 0.
      df(2,3) = -el*cons*z(5)
      df(2,4) = -el*en*2.*z(4) + el**2*s
      df(2,5) = -el*cons*z(3)
      return
      end
c................................................................
      subroutine gsub(i, z, g)
      real z(1), g
      go to (10, 20, 30, 40, 30), i
   10 g = z(1)
      return
   20 g = z(3)
      return
   30 g = z(4)
      return
   40 g = z(1) - 1.
      return
      end
c................................................................
      subroutine dgsub(i, z, dg)
      real z(1), dg(1)
      do 10 j=1,5
        dg(j) = 0.
   10 continue
      go to (20, 30, 40, 20, 40), i
   20 dg(1) = 1.
      return
   30 dg(3) = 1.
      return
   40 dg(4) = 1.
      return
      end