C       PROGRAM NAME:  PROB48.F
C	  MECE 6373 VISCOUS FLOW II      
C       PROGRAM DEVELOPED TO SOLVE PROBLEM 4-8
C       SOLVES UP TO TEN SIMULTANEOUS FIRST ORDER ORDINARY
C       DIFFERENTIAL EQUATIONS USING RUNGE-KUTTA MARCHING SCHEME
C       AS OUTLINED IN F.M. WHITE, VISCOUS FLUID FLOW, 2ND ED, 1991, APP C
C	  DR. CONSTANTINE TARAWNEH (LAST REVISED 2/14/2008)
C
	IMPLICIT NONE
      COMMON/BETAD/ BETA,XLIM     
      REAL Y1, Y10, Y11, Y2, Y20, Y21
      REAL BETA, X, XLIM, H 
	PRINT*,'ENTER BETA,XLIM,Y10'
	READ*,BETA,XLIM,Y10
C	
      CALL TRY(Y10,Y20)
	Y11 = Y10+0.01
C	
	CALL TRY(Y11,Y21)
	Y1 = (Y11-Y10)*(1.0-Y20)/(Y21-Y20)+Y10
 888	PRINT*, ' Y1 = ',Y1	
	Y10 = Y11
	Y11 = Y1
	Y20 = Y21
	CALL TRY(Y10,Y20)
	IF(ABS(Y20-1.00000).LT.0.000001) GOTO 100
	CALL TRY(Y11,Y21)
	Y1 = (Y11-Y10)*(1.0-Y20)/(Y21-Y20)+Y10	
	GOTO 888
 100	PRINT*, ' Y1 = ',Y1
	END
C
	SUBROUTINE TRY(Y1,Y2)
	COMMON/BETAD/ BETA,XLIM
	REAL Y1, Y10, Y11, Y2, Y20, Y21
      REAL BETA, X, XLIM, H 
      DIMENSION Y(10),F(10)
	Y(1)= Y1
	Y(2)= 0.0
	Y(3)= 0.0
	X = 0.0
	N = 3
	H = 0.05
	M = 0
  8	IF(X-XLIM) 6,6,7
  6	CALL RUNGE(N,Y,F,X,H,M,K)
	GOTO (10,20),K
  10	F(1)= -Y(3)*Y(1)-BETA*(1.0-Y(2)*Y(2))
	F(2)= Y(1)
	F(3)= Y(2)
	GOTO 6
  20	CONTINUE
	GOTO 8
  7	Y2 = Y(2)
	RETURN
	END
C
	SUBROUTINE RUNGE(N,Y,F,X,H,M,K)
C	THIS ROUTINE PERFORMS RUNGE-KUTTA CALCULATION BY 
C	GILLS METHOD
C	F.M. WHITE, VICOUS FLUID FLOW, 2ND ED, 1991,APP C
C      
      REAL Y1, Y10, Y11, Y2, Y20, Y21
      REAL BETA, X, XLIM, H 
	DIMENSION Y(10), F(10), Q(10)
	M=M+1
	GOTO (1,4,5,3,7),M
  1	DO 2 I=1,N
  2	Q(I) =0.0
	A=0.5
	GOTO 9
  3	A=1.707107
C	IF YOU NEED MORE ACCURACY, USE
C	A=1.7071067811865475244
  4	X=X+0.5*H
  5	DO 6 I=1,N
	Y(I)=Y(I)+A*(F(I)*H-Q(I))
  6	Q(I)=2.0*A*H*F(I)+(1.0-3.0*A)*Q(I)
	A=0.2928932
C	IF YOU NEED MORE ACCURACY, SET
C	A=0.2928932188134524756
	GOTO 9
  7	DO 8 I=1,N
  8	Y(I)=Y(I)+H*F(I)/6.0-Q(I)/3.0
	M=0
	K=2
	GOTO 10
  9	K=1
  10	RETURN
	END