! Driver program for SNK14-Ex2  
! The BVODE to be solved, on the interval [0,1], is 
! y_1' = y_2, y_2' = (-2*x^2 - x^3/6 - 1/(4+x)^2) + x*y_3 - y_4, y_3' = exp(y_1), y_4' = x*exp(y_1)
! y_1(0)=ln(4),  y_3(0) = 0, y_4(0) = 0,  y_1(1)=ln(5).
! Exact solution is y_1(x) = ln(4+x)

       MODULE DEFINE_FCN

!---------------------------------------------------------------------------------------

        IMPLICIT NONE

        ! Global variables and constants. 


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

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

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


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

        ! The user supplied subroutines FSUB, BCSUB, and GUESS. 

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

         ! Input arguments
         DOUBLE PRECISION::X ! point where right hand side function is evaluated  
         DOUBLE PRECISION::Z(NSOL) ! solution approximation vector at X
 
         ! Output arguments
         DOUBLE PRECISION::F(NODE) ! function values

         F(1)= Z(2)
         F(2)= (-2.0D0*X**2-X**3/6.0d0-1.0d0/(4.0d0+X)**2)+X*Z(3)-Z(4)
         F(3)= exp(Z(1))
         F(4)= X*exp(Z(1))
         
       END SUBROUTINE FSUB



       ! Subroutine to evaluate the boundary conditions

       SUBROUTINE BCSUB(I,Z,G)

        ! Input Arguments
        INTEGER :: I ! Ith boundary condition, starting with the left boundary
        DOUBLE PRECISION :: Z(NSOL) ! solution approximation vector at boundary point

        ! Output Arguments
        DOUBLE PRECISION :: G ! Boundary values


        IF (I.EQ.1) THEN
          G = Z(1) - log(4.0d0)
        END IF

        IF (I.EQ.2) THEN
          G = Z(3)
        END IF

       IF (I.EQ.3) THEN
          G = Z(4)
        END IF

        IF (I.EQ.4) THEN
          G = Z(1)-log(5.0d0)
        END IF

 
       END SUBROUTINE BCSUB


       SUBROUTINE GUESS(X, Z, DMVAL)
  
        ! Input arguments
        DOUBLE PRECISION::X ! Point where solution guess is required

        ! Output Arguments
        DOUBLE PRECISION:: Z(NSOL) ! Initial approximation values
        DOUBLE PRECISION:: DMVAL(NODE)  ! Initial approximation derivative values

        ! Local variables
        INTEGER :: I  ! Loop counter

     
         DO I = 1,NODE
          DMVAL(I) = 0.D0
         ENDDO

         Z(1) = X+1.0d0
         Z(2) = 1.0d0
         Z(3) = X**2/2.0D0 + x
         Z(4) = X*Z(3)


        END SUBROUTINE GUESS



       END MODULE DEFINE_FCN

       !----------------------------------------------------------------------------------------------------
       ! Main program 
       !----------------------------------------------------------------------------------------------------
       PROGRAM SNK14Ex2

       USE DEFINE_FCN

       USE COLNEW_WRAPPER

       IMPLICIT NONE

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

       INTEGER :: I,J ! Loop counters
       DOUBLE PRECISION :: T, Z(NSOL)  ! Evaluation point + approximate solution
       DOUBLE PRECISION:: X(2) ! Initial mesh; uniform mesh of 5 subinterval will be
                                ! automatically generated
       DOUBLE PRECISION :: A= 0.0D0, B=1.0D0 ! Left and right endpoints  
       DOUBLE PRECISION :: TOL ! User tolerance 
       INTEGER :: K_POINTS ! Number of collocation points

       ! Set FORMAT
       CHARACTER(LEN=*), PARAMETER  :: FMT = "(F3.1,A3,ES10.3,A3,ES10.3,A3,ES10.3,A3)"
      
       !---------------------------------------------------------------------------------------------------------
      
     
        ! Set initial mesh array; uniform mesh of 5 subinterval will be automatically generated
        X=(/A,B/) 

        ! Set the value for the number of collocation points
        K_POINTS=4

        ! Set the tolerance.       
        TOL = 1.0D-6


        ! Initialize SOL data structure
        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)                         
          

         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)

           ! Print out all solution components
           DO J=1,NSOL
            PRINT*,'THE SOLUTION APPROXIMATION OF COMPONENT ',J, &
                    'FOR T= ',T,'IS Z = ',Z(J)
           ENDDO
           !WRITE(*,*)

           ! Print out first solution approximation and its error
           !PRINT*,'THE SOLUTION APPROXIMATION OF COMPONENT ',1, &
           !         'FOR T= ',T,'IS Z = ',Z(1),' ERROR = ', (ABS(Z(1)-log(4.0d0+T)))
           !WRITE(*,*)

           ! Print out info for table
           !WRITE(*,FMT) T, " & ", log(4.0d0+T), " & ", Z(1), " & ", &
	   !				ABS(Z(1)-log(4.0d0+T)), " \\"
 
          ENDDO

         endif
          
         ! The function colnew_terminate is used to free up dynamic memory            
         CALL COLNEW_TERMINATE(SOL)    
       

       END PROGRAM SNK14Ex2


