c     Problem 3 for colnew, example of continuation 


c     Quantum Neumann equation, with 2 "eigenvalues" c_1 and c_2     
c     (c_1=v-c_2-c_3, v is a parameter, used in continuation)        
c                                                                
c     diff(f,x,2) + (1/2)*(1/x + 1/(x-1) + 1/(x-y))*diff(f,x)        
c      - (c_1/x + c_2/(x-1) + c_3/(x-y))* f(x) = 0               
c                                                                
c     with 4 "boundary" conditions: diff(f,x)(a_k)=2*c_k*f(a_k) for   
c     k=1,2,3, a_k=(0, 1, y)  and normalization f(1) = 1             
c                                                                
c     The z-vector is z_1=f, z_2=diff(f,x), z_3=c_2 and z_4=c_3      
c     The guess is chosen to have one node in [0,1],  f(x)=2*x-1     
c     such that f(1)=1, c_2 and c_3 are chosen to cancel poles in    
c     the differential equation at 1.0 and y, z_3=1, z_4=1/(2*y-1)   


      subroutine fsub(x, z, f)
      common /cprob3/v, y
      real*8 v, y
      real*8 z(4), f(3), x      
      f(1) = - (1/2.0d0)*(1/x + 1/(x-1.0d0) + 1/(x-y))*z(2) +
     $     ((v - z(3) - z(4))/x + z(3)/(x-1.0d0) + z(4)/(x-y))*z(1)
      f(2) = 0
      f(3) = 0
      return
      end
      subroutine dfsub(x, z, df)
      common /cprob3/v, y
      real*8 v, y
      real*8 z(4), f(3), x, df(3, 4)
      df(1,1) = z(4)/(x-y)+(-z(4)-z(3)+v)/x+z(3)/(x-1.0d0)
      df(1,2) =-0.5d0*(1/x+1/(x-1.0d0)
     $     +1/(x-1.9d0))
      df(1,3) = (1/(x-1.0d0)-1/x)*z(1)
      df(1,4) =(1/(x-y)-1/x)*z(1)
      df(2,1) =0
      df(2,2) =0
      df(2,3) =0
      df(2,4) =0
      df(3,1) =0
      df(3,2) =0
      df(3,3) =0
      df(3,4) =0
c     print *, '  df =', df
      return
      end
      subroutine gsub(i, z, g)
      common /cprob3/v, y
      real*8 v, y
      real*8 z(4), g
      go to (10, 20, 30, 40), i
 10   g = z(2) - 2.0d0*z(1)*(v-z(3)-z(4))
      go to 99
 20   g = z(2) - 2.0d0*z(1)*z(3)
      go to 99
 30   g = z(1) - 1.0d0
      go to 99
 40   g = z(2) - 2.0d0*z(1)*z(4)
 99   continue
      return
      end
      subroutine dgsub(i, z, dg)
      common /cprob3/v, y
      real*8 v, y
      real*8 z(4), dg(4)
      go to (10, 20, 30, 40), i
 10   dg(1) = -2.0d0*(-z(4)-z(3)+v)
      dg(2) = 1.0d0
      dg(3) = 2.0d0*z(1)
      dg(4) = 2.0d0*z(1)
      go to 99
 20   dg(1) = -2.0d0*z(3)
      dg(2) = 1.0d0
      dg(3) = -2.0d0*z(1)
      dg(4) = 0
      go to 99
 30   dg(1) = 1.0d0
      dg(2) = 0
      dg(3) = 0
      dg(4) = 0
      go to 99
 40   dg(1) = -2.0d0*z(4)
      dg(2) = 1.0d0
      dg(3) = 0
      dg(4) = -2.0d0*z(1)
 99   continue
      return
      end

      subroutine guess(x, z, dmval)
      common /cprob3/v, y
      real*8 v, y
      real*8 x, z(4), dmval(3)
      z(1) = 2.0d0*x-1.0d0
      z(2) = 2.0d0
      z(3) = 1.0d0
      z(4) = 1/(2.0d0*y-1.0d0)
      dmval(1) = 0
      dmval(2) = 0
      dmval(3) = 0
      return
      end
      
      program prob3
      implicit real*8 (a-h, o-z)
      real*8 fspace(30000), zeta(4), tol(4), z(4)
      integer ispace(2000), m(3), ipar(11), ltol(4)            
      external fsub, dfsub, gsub, dgsub, guess                       

      real*8 fixpnt(5)
      common /cprob3/v, y
      real*8 v, y

      parameter (nsteps=100)
      real*8 eig(nsteps, 4)
      
      y = 1.9d0
      v = 0.00001d0
      m(1) = 2
      m(2) = 1
      m(3) = 1

      zeta(1) = 0.0d0
      zeta(2) = 1.0d0
      zeta(3) = 1.0d0
      zeta(4) = y

      do 10 i = 1, 11
         ipar(i) = 0
 10   continue

c     Non-linear prob
      ipar(1) = 1
c     Number of collocation points
      ipar(2) = 3
c     Initial uniform mesh of 4 subintervals
      ipar(3) = 4
      ipar(8) = 0
c     Size of fspace, ispace
      ipar(5) = 30000
      ipar(6) = 2000
c     medium output, set it to 1 to have a clean table of eigenvalues
      ipar(7) = 0
c     Initial approx provided
      ipar(9) = 1
c     fixpnt is an array containing all fixed points in the mesh
      ipar(11) = 1
c     tolerances on all components
      ipar(4) = 4
c     tolerance on f and diff(f,x) and on c2 and c3
      do 20 i = 1, 4
         ltol(i) = i
         tol(i) = 1d-5
 20   continue
c     only one interior boundary
      fixpnt(1) = 1.0d0

c     First run, using the guess
      call colnew(3, m, 0d0, y, zeta, ipar, ltol, tol, fixpnt, ispace,
     $     fspace, iflag, fsub, dfsub, gsub, dgsub, guess)
      len = int(1+y/.1d0)
      print *, 'Wave function at v=0'
      do 30 i = 0, len
         call appsln(.1d0*i, z, fspace, ispace)
c     Print z(1) to see the solution here, a straight line
      print *, z(1)
 30   continue
      print *, 'Eigenvalues'
      write(*, 9999) v, z(3), z(4), z(3)+y*z(4)

c     do continuation, which implies modifying ipar(9) and ipar(3)
      ipar(9) = 3
      ipar(8) = 1
      do 100 i = 1, nsteps
         ipar(3) = ispace(1)
c     poor man's logarithmic scale
         if ( i .lt. nsteps/3) v = v + 0.2d0
         if ( i .ge. nsteps/3) v = v + 0.4d0
         if ( i .gt. 2*nsteps/3) v = v + 0.6d0
         call colnew(3, m, 0d0, y, zeta, ipar, ltol, tol, fixpnt,
     $        ispace,
     $        fspace, iflag, fsub, dfsub, gsub, dgsub, guess)
         call appsln(.1d0, z, fspace, ispace)
         eig(i,1) = z(3)
         eig(i,2) = z(4)
         eig(i,3) = z(3)+y*z(4)
         eig(i,4) = v
 100  continue

      do 110 i = 1, nsteps
         write(*,9999) eig(i, 4), eig(i, 1), eig(i, 2), eig(i, 3)
 110  continue
 9999 format(5H v = ,f6.3,7H c_2 = ,f6.3,7H c_3 = ,f6.3,7H E/2 = ,f6.3)

c     Show the solution for the last value of v
      print *, 'Wave function at large v'
      do 120 i = 0, len
         call appsln(.1d0*i, z, fspace, ispace)
         print *, z(1)
 120  continue

c     The results show that the probability of presence is concentrated
c     near t=1 and near t=y which corresponds to a pole of the sphere.

      
      end