C
C PAGE 421-420: NUMERICAL MATHEMATICS AND COMPUTING, CHENEY/KINCAID, 1985
C
C FILE: BVP2.FOR
C
C BOUNDARY VALUE PROBLEM SOLVED BY SHOOTING METHOD (RK4SYS,XPSYS)
C
      DIMENSION  X(5),X2(99),X4(99)   
      DATA  NP1/99/, TA/1.0/, TB/2.0/ 
      DATA  ALPHA/1.09737491/, BETA/8.63749661/ 
      DATA  X/1.0,1.09737491,0.0,1.09737491,1.0/
      H = (TB - TA)/REAL(NP1) 
      T = TA      
      DO 2 I = 1,NP1
        CALL RK4SYS(5,T,X,H,1,0)      
        X2(I) = X(2)
        X4(I) = X(4)
        T = TA + REAL(I)*H  
   2  CONTINUE
      P = (BETA - X4(NP1))/(X2(NP1) - X4(NP1))  
      Q = 1.0 - P 
      DO 3 I = 1,NP1
        X2(I) = P*X2(I) + Q*X4(I)     
   3  CONTINUE
      ERROR = EXP(TA) - 3.0*COS(TA) - ALPHA     
      PRINT 5,TA,ALPHA,ERROR
      DO 4 I = 9,NP1,9      
        T = TA + REAL(I)*H  
        ERROR = EXP(T) - 3.0*COS(T) - X2(I)     
        PRINT 5,T,X2(I),ERROR 
   4  CONTINUE
   5  FORMAT(5X,F10.5,2(5X,E22.14))   
      STOP
      END 
        
      SUBROUTINE XPSYS(X,F) 
      DIMENSION  X(5),F(5)  
      F(1) = 1.0  
      F(2) = X(3) 
      F(3) = EXP(X(1)) - 3.0*SIN(X(1)) + X(3) - X(2)      
      F(4) = X(5) 
      F(5) = EXP(X(1)) - 3.0*SIN(X(1)) + X(5) - X(4)      
      RETURN      
      END 
        
      SUBROUTINE RK4SYS(N,T,X,H,NSTEP,J)
      DIMENSION  X(10),Y(10),F1(10),F2(10),F3(10),F4(10)  
      H2 = 0.5*H  
      START = T   
      DO 6 K = 1,NSTEP      
      CALL XPSYS(X,F1)      
      DO 2 I = 1,N
   2  Y(I) = X(I) + H2*F1(I)
      CALL XPSYS(Y,F2)      
      DO 3 I = 1,N
   3  Y(I) = X(I) + H2*F2(I)
      CALL XPSYS(Y,F3)      
      DO 4 I = 1,N
   4  Y(I) = X(I) + H*F3(I) 
      CALL XPSYS(Y,F4)      
      DO 5 I = 1,N
   5  X(I) = X(I) + H*(F1(I) + 2.0*(F2(I) + F3(I)) + F4(I))/6.0     
      T = START + REAL(K)*H 
   6  IF(K .EQ. (K/J)*J)  PRINT 7,T,(X(I),I = 1,N)
   7  FORMAT(5X,F10.5,4(5X,E22.14))   
      RETURN      
      END