        ! Driver program for TP20 from the bvp test website
        ! This problem uses the analytic derivatives, DFDY and DBCDY
        


       MODULE DEFINE_FCN


        IMPLICIT NONE
           
        DOUBLE PRECISION::EPS=0.05d0
        ! Problem defining parameter

        INTEGER,PARAMETER::NODE=2
        ! Node is the number of differential equations

        INTEGER,PARAMETER::NSOL=2
        ! Nsol is the number of solution and derivative components. 

        INTEGER,PARAMETER::LEFTBC=1
        ! Leftbc is the number of boundary conditions at the left side boundary point.

        INTEGER,PARAMETER,DIMENSION(2)::ORDERS=(/1,1/)
        ! Orders is an array containing the order of each differential equation
       
       CONTAINS


         ! Subroutine to evaluate the function
        SUBROUTINE FSUB(X,Z,F)

         ! Input arguments
         DOUBLE PRECISION::X ! meshpoint 
         DOUBLE PRECISION::Z(NSOL) ! solution approximation

         ! Output arguments
         DOUBLE PRECISION::F(NODE) ! function values


         F(1)=Z(2)
         F(2)= (1.0D0 - (Z(2)**2) ) / EPS
    
   
        END SUBROUTINE FSUB


        SUBROUTINE DFDY(X,Z,DF)

         ! Input arguments
         DOUBLE PRECISION :: X, Z(NSOL)

         ! Output Argument
         DOUBLE PRECISION :: DF(NODE,NSOL)

         ! Local Variables
         INTEGER :: I, J
       
         DO I = 1,NODE
           DO J = 1,NSOL
            DF(I,J) = 0.D0
           ENDDO
         ENDDO

         DF(1,1)=0.0D0
         DF(1,2)=1.0D0
       
         DF(2,1) = 0.0d0
         DF(2,2) = (-2.0d0*Z(2))/EPS
       
       END SUBROUTINE DFDY

       ! Subroutine to evaluate the boundary conditions

       SUBROUTINE BCSUB(I,Z,G)

         ! Input Arguments
         INTEGER :: I ! Ith boundary condition, starting with the left end boundary
         DOUBLE PRECISION :: Z(NSOL)

         ! Output Arguments
         DOUBLE PRECISION :: G ! Boundary values

         IF (I.EQ.1) THEN
          G = Z(1) - 1.0D0 - EPS*  &
                    LOG(COSH(-0.745D0/EPS))
         END IF
         IF (I.EQ.2) THEN 
          G = Z(1) - 1.0D0 - EPS*  &
                    LOG(COSH(0.255D0/EPS))
         END IF

      
       END SUBROUTINE BCSUB
      

      SUBROUTINE DBCDY(I,Z,DG)

        !Input Arguments
        INTEGER :: I
        DOUBLE PRECISION :: Z(NSOL)

        !Output Arguments
        DOUBLE PRECISION :: DG(NSOL)

        !Local Variables
        INTEGER :: J
        
        DO J = 1,NSOL
            DG(J) = 0.0D0
        ENDDO
        IF ((I.EQ.1).OR.(I.EQ.2)) THEN
            DG(1) = 1.D0
        END IF
        
       END SUBROUTINE DBCDY

       SUBROUTINE GUESS(X, Z, DMVAL)

        !Input Arguments
        DOUBLE PRECISION :: X

        !Output Arguments
        DOUBLE PRECISION :: Z(NSOL), DMVAL(NODE)
        !Local Variables
        INTEGER :: I
       
        DO I = 1,NODE
         DMVAL(I) = 0.D0
       ENDDO

       Z(1) = 1.0D0
       Z(2) = 1.0D0

      END SUBROUTINE GUESS

     END MODULE DEFINE_FCN

       !----------------------------------------------------------------------------------------------------
       ! Start of the program which drives all functions and subroutines in the two modules define_fcn and
       ! colnew_wrapper
       !----------------------------------------------------------------------------------------------------
       PROGRAM CASH20_FIRSTORDER

        USE DEFINE_FCN

        USE COLNEW_WRAPPER

        IMPLICIT NONE

        TYPE(COLNEW_SOL)::SOL ! Stores most of the information specific to the final solution.
      

       INTEGER :: I,J
       DOUBLE PRECISION :: T, Z(NSOL) 
       DOUBLE PRECISION:: X(2)! Initial mesh
       DOUBLE PRECISION :: A= 0.0D0, B=1.0D0 ! Left and right endpoints  
       DOUBLE PRECISION :: TOL ! User tolerance 
       INTEGER :: K_POINTS !  Define the number of collocation points
      
       !---------------------------------------------------------------------------------------------------------
      
     

        ! A uniform mesh will be constructed, starting at the left endpoint and ending at the right endpoint

        X=(/A,B/) 

        ! Setting the value for the number of collocation points

        K_POINTS=7

        ! Setting the value for the user tolerance.       

        TOL = 1.0D-8


        ! Both the hidden colnew_init and colnew_solver functions are necessary to compute a solution approximation.
        !
        ! The call to colnew_init has to be of the form:
        !
        ! sol = colnew_init(node,leftbc,orders,x,nfixpts)
        !
        ! max_num_intervals, the maximum number of subintervals allowed to compute a solution approximation,  is an optional integer parameter
            
        

        SOL=COLNEW_INIT(NODE,LEFTBC,ORDERS,X,&   
                          MAX_NUM_SUBINTERVALS=10000)


        ! The call to colnew_solver has to be of the form:
        ! 
        ! sol=colnew_solver(sol,fsub,bcsub).
        !
        ! The optional parameters are :
        !
        ! a) subroutine dfdy, the derivative information for fsub
        !
        ! b) subroutine dbcdy, the derivative information for gsub
        !
        ! c) subroutine guess, the initial guess for the solution approximation
        !
        ! The default for the remaining optional parameters are listed below
        ! 
        !
        ! d) nonlinear (integer) = 0  if the problem is linear 
        !              = 1  if the problem is nonlinear (default)
        !
        ! e) num_coll_points (integer) = number of collocation points per subinterval (=k)
        !                      where max(orders).le.k.le 7. 
        !
        !                      k = 3 (default)
        !
        ! f) trace (integer) = output control
        !          = -1 for full diagnostic printout (default)
        !          = 0 for selected printout
        !          = 1 for no printout
        !
        ! g) problem_type (integer) = 0 if the problem is regular (default)
        !                 = 1 if the first relax factor is =rstart, and the nonlinear
        !                     iteration does not rely on past convergence (use for only
        !                     extra sensitive nonlinear problems)
        !                 = 2 if we are to return immediately upon a) two successive
        !                     nonconvergences, or b) after obtaining an error estimate 
        !                     for the first time.
        !
        ! 5) fixpoints (integer) = an array of size "nfixpts" containing the points, other than
        !                the left and right boundary points, which are to be included 
        !                in every mesh.
        !                fixpoints is an empty array by default
        !
        ! 6) ttol (double precision) = tolerance applied to all components of the solution vector.
        !           ttol=1.0d-5 is the default.
        !
        ! 7) mesh_halving (integer) = 1 if initial mesh is provided by the user and adaptive
        !                   mesh selection is to be done (default)
        !                 = 2 if initial mesh is provided by the user but in addition
        !                     NO adaptive mesh selection is to be done.
               
       
         SOL=COLNEW_SOLVER(SOL,FSUB,BCSUB,TTOL=TOL,&
                         NUM_COLL_POINTS=K_POINTS, &
                         GUESS=GUESS)                         
          

        ! The function colnew_terminate is used to free up dynamic memory




         if(sol%info.eq.1)then


          PRINT*,' '         
          PRINT*,'SAMPLE OUTPUT'
          PRINT*,'  '

          DO I = 1,11
           T=A+((I-1)/10.0D0)*(B-A)
           CALL COLNEW_EVAL(SOL,T,Z,FSUB)
           DO J=1,NSOL
            PRINT*,'THE SOLUTION APPROXIMATION OF COMPONENT ',J, &
                    'FOR T= ',T,'IS Z = ',Z(J)
           ENDDO
          ENDDO

         endif
                       
         CALL COLNEW_TERMINATE(SOL)    

     

  
       END PROGRAM CASH20_FIRSTORDER


