C
C ********************************************************************
C
C     Two Point Boundary Value Problem due to B.A. Troesch
C
C     y''(t) = xlam * sinh(xlam * y(t))
C     y(0) = 0, y(1) = 1
C
C     J.Stoer, R.Bulirsch: Introduction to Numerical
C     Analysis, Springer, 2002 (Third Edition)
C
C     Here: xlam = 5
C
C ********************************************************************
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      PARAMETER (MMAX=21,MSMAX=2,MMS=23,NMS=6,
     *           NDW=NMS*(8+MMAX*(10+6*NMS)),NDIW=5*NMS)
      PARAMETER (NFILE=6)
      DIMENSION X(MMS),XS(MSMAX),Y(NMS,MMS),WORK(NDW)
      INTEGER JS(MMS,MSMAX),IWORK(NDIW)
      EXTERNAL F, R, DIFSYB, RKF7B
      COMMON /ITEIL/ ITEIL
      COMMON /PARAM/ XLAM
C
      OPEN(6,FILE='troesch.txt',STATUS='UNKNOWN')
      OPEN(7,FILE='troesch.dat',STATUS='UNKNOWN')
C
C --- Parameters
C
      N     =  2
      M     =  21
      MS    =  0
      TOL   =  1.D-08
      ITMAX =  20
      KS    =  0
C
C --- Initial Trajektory
C
      DO 100 J=1,M
      X(J) = DFLOAT(J-1)/DFLOAT(M-1)
      Y(1,J) = 0.D0
      Y(2,J) = 0.D0
  100 CONTINUE
C
C --- Parameter lambda of the  Troesch Problem
C     (may be changed)
C
      XLAM  = 5.0D0
      WRITE(6,1000) XLAM
C
C --- Solution of the Boundary Value Problem
C
      ITEIL = 1
      KP    = 0
      CALL BNDSCO(F,R,DIFSYB,X,XS,Y,WORK,IWORK,JS,N,M,MS,
     *     KS,TOL,ITMAX,KP,MMAX,MSMAX,MMS,NMS,NDW,NDIW,NFILE)
      IF(KP.LT.0) GOTO 900
C
C --- Numerical Integration of the Solution Trajectory,
C     Determination of the Data for Plot with MATLAB
C     Files: troesch.dat, troesch.m
C 
      ITEIL  =  2
      KP     =  0
      IFEIN  =  5
      NPLOT  =  7
      IPLOT  =  1
      ICASE  =  1
      CALL AWP(DIFSYB,F,X,XS,Y,N,M,MS,
     *         JS,MMS,NMS,IFEIN,IPLOT,KP,NFILE,NPLOT,ICASE)
C
  900 CONTINUE
 1000 FORMAT(///' BVP of Troesch: Lambda =',D15.5)
C -------------------------------------------------------------------
      STOP
      END



      SUBROUTINE F(X,Y,DY,J,L,JS,MMS)
C
C     Two Point Boundary Value Problem due to B.A. Troesch
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DIMENSION Y(*),DY(*)
      INTEGER JS(MMS,*)
      COMMON /ITEIL/ ITEIL
      COMMON /PARAM/ XLAM
C
      DY(1) = Y(2)
      DY(2) = XLAM*SINH(XLAM*Y(1))
C
C -------------------------------------------------------------------
      RETURN
      END




      SUBROUTINE R(YA,YB,ZZ,W,NYA,NSK,J,L,LS,JS,MMS)
C
C     Two Point Boundary Value Problem due to B.A. Troesch
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DIMENSION YA(*),YB(*),ZZ(*),W(*)
      INTEGER NYA(*),NSK(*)
      INTEGER JS(MMS,*)
      COMMON /ITEIL/ ITEIL
      COMMON /PARAM/ XLAM
C
      W(1)   =  YA(1)
      W(2)   =  YB(1)-1.D0
C
      NYA(1) =  1
C -------------------------------------------------------------------
      RETURN
      END




      SUBROUTINE AWP(METHOD,FCN,X,XS,Y,N,M,MS,JS,MMS,NMS,IFEIN,
     *               IPLOT,KP,NFILE,NPLOT,ICASE)
C
C *******************************************************************
C
C     INTEGRATION OF INITIAL VALUE PROBLEMS
C
C     IFEIN : NUMBER OF INTEGRATION NODES IN SUBINTERVAL
C     IPLOT : THE NODES ARE STORED FOR PLOT SUBROUTINE,
C             IF IPLOT.EQ.1 THE ROUTINE TAKES IFEIN = 40
C             EQUIDISTANT POINTS IN EACH SUBINTERVAL
C     KP    : PRINT PARAMETER
C             KP.LT.0 : NO PRINT
C             KP.EQ.0 : MULTIPLE SHOOTING NODES
C             KP.GT.0 : ALL INTEGRATION NODES
C     ICASE : 1  :  ROUTINE INTEGRATES OUTPUT OF BNDSCO IN EACH
C                   SUBINTERVAL (THE SWITCHING POINTS ARE ALREADY
C                   ARRANGED IN ORDER; THE ROUTINE TAKES NEW INITIAL
C                   DATA FOR EACH SUBINTERVAL)
C             0  :  INTEGRATION OVER THE WHOLE INTERVAL; THE ALGORITHM
C                   ARRANGES THE SWITCHING POINTS XS(K) IN ORDER AND
C                   SAVES THE DATA Y(I;J);
C
C *******************************************************************
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DIMENSION X(*),XS(*),Y(NMS,*)
      DIMENSION Z(20),DZ(20)
      INTEGER JS(MMS,*)
      EXTERNAL FCN , R
      COMMON /ITEIL/ ITEIL
      COMMON /PARAM/ XLAM
C
      IF(IPLOT.GT.0) IFEIN=40
      TOL =  1.D-07
C
      NN  =  N
C
C     ARRANGING THE SWITCHING POINTS (IN CASE ICASE = 0)
C
      IF(MS.EQ.0.OR.ICASE.GT.0) GOTO 300
      K   =  1
  200 CONTINUE
      K1  =  K + 1
      X1  =  X(K)
      X2  =  X(K1)
      DO 210 L=1,MS
      IF((X1-XS(L))*(X2-XS(L)).LT.-1.D-10) GOTO 230
  210 CONTINUE
      GOTO 250
  230 CONTINUE
      M   =  M + 1
      DO 240 J=M,K+2,-1
      X(J) = X(J-1)
  240 CONTINUE
      X(K1) = XS(L)
      GOTO 200
  250 K   =  K1
      IF(K.LT.M) GOTO 200
  300 CONTINUE
C
      JT  =  0
      L   =  0
      HI  =  0.D0
      M1  =  M - 1
      IF(KP.GE.0) THEN
        WRITE(NFILE,2000)
        WRITE(NFILE,2200)
      ENDIF
C
C     LOOP OVER THE SUBINTERVALS
C
      DO 500 K=1,M1
      J=K
      J1=J+1
      L=L+1
C
      T=X(J)
      TEND=X(J1)
      IF(MS.GT.0) THEN
         DO 503 KL=1,MS
  503    IF(T.EQ.XS(KL)) L=L-1
      ENDIF
      HMAX=(TEND-T)/IFEIN
      IF(HI.EQ.0.D0) HI=0.01D0*HMAX
      HMAX=ABS(HMAX)
C
C     DEFINITION OF JS(J,K), K=1,..,MS (ICASE = 0)
C
      IF(ICASE.EQ.0) THEN
      DO 410 KL=1,MS
      JS(J,KL) = 1
      IF(T.LT.XS(KL)-1.D-10) JS(J,KL) = -1
  410 CONTINUE
      END IF
C
C     NEW INITIAL DATA FOR THE INTEGRATION
C
      IF(ICASE.EQ.0.AND.J.GT.1) GOTO 430
      DO 420 I=1,N
  420 Z(I)=Y(I,J)
  430 CONTINUE
C
      CALL FCN(T,Z,DZ,J,L,JS,MMS)
C
      IF(KP.GE.0) WRITE(NFILE,2001)
      JT=JT+1
      IF(KP.GE.0)
     *WRITE(NFILE,3000) T,(Z(I),I=1,NN)
      IF(IPLOT.EQ.1)
     *WRITE(NPLOT,4000) T,(Z(I),I=1,NN)
C
C     LOOP FOR SUB-SUB-INTERVALS
C
      DO 510 K1=1,IFEIN
      T1=T+SIGN(HMAX,HI)
      IF(K1.EQ.IFEIN) T1=TEND
C
      CALL METHOD(NN,FCN,T,Z,T1,TOL,HMAX,HI,J,L,JS,MMS)
C
      IF(HI.EQ.0.D0) GOTO 600
C
      CALL FCN(T,Z,DZ,J,L,JS,MMS)
C
      JT=JT+1
      IF(KP.GT.0)
     *WRITE(NFILE,3000) T,(Z(I),I=1,NN)
      IF(IPLOT.EQ.1)
     *WRITE(NPLOT,4000) T,(Z(I),I=1,NN)
C
  510 CONTINUE
C
C     DEFINITION OF THE MULTIPLE SHOOTING KNOTS (ICASE = 0)
C
      IF(ICASE.EQ.0) THEN
      DO 520 I=1,NN
      Y(I,J+1) = Z(I)
  520 CONTINUE
      ENDIF
C
  500 CONTINUE
C
      IF(KP.EQ.0) THEN
      WRITE(NFILE,2001)
      WRITE(NFILE,3000) T,(Z(I),I=1,NN)
      ENDIF
      IF(KP.GE.0) WRITE(NFILE,2001)
      GOTO 900
C
  600 CONTINUE
      WRITE(NFILE,5000)
      WRITE(NFILE,3000) T,(Z(I),I=1,NN)
C
  900 CONTINUE
 1000 FORMAT(I4)
 2000 FORMAT(///)
 2001 FORMAT(' ',77('*'))
 2200 FORMAT(7X,'t',17X,'y',17X,'dy/dt')
 3000 FORMAT(4D18.10/5(18X,3D18.10/))
 4000 FORMAT(5(30E13.5/))
 5000 FORMAT(' Termination: Integration stepsize H = 0')
C -------------------------------------------------------------------
      RETURN
      END



      SUBROUTINE RESET(X,XS,Y,N,M,MS,NMS)
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DIMENSION X(*),XS(*),Y(NMS,*)
C
      J1=0
      MD=M
      DO 600 J2=1,M
      DO 601 I=1,MS
      IF(XS(I).EQ.X(J2)) GOTO 620
  601 CONTINUE
      J1=J1+1
      X(J1)=X(J2)
      DO 610 I=1,N
  610 Y(I,J1)=Y(I,J2)
      GOTO 600
  620 MD=MD-1
  600 CONTINUE
      M=MD
C -------------------------------------------------------------------
      RETURN
      END



