
      SUBROUTINE INITAL
C...
C...  M1 - BLASIUS PROBLEM
C...
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      COMMON/T/    ETA,  NSTOP,  NORUN
     +      /Y/     G1,     G2,     G3
     +      /F/    G1E,    G2E,    G3E
     +      /S/    G3S,   DG3S,  FLAGM,  FLAGP
     +      /I/     IP
C...
C...  SET THE INITIAL ESTIMATE OF G3, INCREMENT IN G3,
C...  CONTROL VARIABLES FOR ITERATION ON G3
      IF(NORUN.EQ.1)THEN
         G3S =0.3D0
         DG3S=0.2D0
         FLAGM=-1.0
         FLAGP=-1.0
      END IF
C...
C...  INITIAL CONDITIONS
      G1=0.0D0
      G2=0.0D0
      G3=G3S
C...
C...  INITIAL DERIVATIVES
      CALL DERV
      IP=0
      RETURN
      END
      SUBROUTINE DERV
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      COMMON/T/    ETA,  NSTOP,  NORUN
     +      /Y/     G1,     G2,     G3
     +      /F/    G1E,    G2E,    G3E
     +      /S/    G3S,   DG3S,  FLAGM,  FLAGP
     +      /I/     IP
C...
C...  ODES
      G1E=G2
      G2E=G3
      G3E=-0.5D0*G1*G3
      RETURN
      END
      SUBROUTINE PRINT(NI,NO)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      COMMON/T/    ETA,  NSTOP,  NORUN
     +      /Y/     G1,     G2,     G3
     +      /F/    G1E,    G2E,    G3E
     +      /S/    G3S,   DG3S,  FLAGM,  FLAGP
     +      /I/     IP
C...
C...  INCREMENT THE COUNTER TO DETERMINE THE END OF EACH RUN
      IP=IP+1
C...
C...  PRINT THE SOLUTION FOR THE FIRST AND LAST RUNS
      IF((NORUN.EQ.1).OR.(NORUN.EQ.20))THEN
C...
C...     PRINT A HEADING
         IF(IP.EQ.1)WRITE(NO,1)NORUN
1        FORMAT(//,' NORUN = ',I3,//,
     +          7X,'eta',8X,'g1',8X,'g2',8X,'g3')
C...
C...     PRINT THE SOLUTION
         WRITE(NO,2)ETA,G1,G2,G3
2        FORMAT(F10.2,3F10.5)
      END IF
C...
C...  AT THE END OF EACH RUN, INCREMENT G3 DEPENDING ON THE SIGN
C...  DEVIATION FROM THE PRESCRIBED BOUNDARY CONDITION G2(INF)
C...  = 1
      IF(IP.EQ.11)THEN
         CALL SEARCH
C...
C...     MONITOR THE ITERATION
         WRITE(*,3)NORUN,G2,G3S,DG3S,FLAGP,FLAGM
3        FORMAT(I5,3F12.6,2F10.3)
      END IF
      RETURN
      END
      SUBROUTINE SEARCH
C...
C...  SUBROUTINE SEARCH SEARCHES FOR A ROOT OF A FUNCTION BY
C...  INTERVAL HALVING
C...
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      COMMON/T/    ETA,  NSTOP,  NORUN
     +      /Y/     G1,     G2,     G3
     +      /F/    G1E,    G2E,    G3E
     +      /S/    G3S,   DG3S,  FLAGM,  FLAGP
     +      /I/     IP
C...
C...  IF (1 - G2(INF)) IS NONNEGATIVE, INCREASE G3(0)
      IF(FUNC(G2).GE.0.0D0)THEN
C...
C...     IF A SMALLER INCREMENT IN G3(0) IS REQUIRED, HALVE
C...     DG3S
         IF(FLAGM.GE.0.0D0)DG3S=DG3S/2.0D0
C...
C...     INCREMENT G3(0) POSITIVELY FOR THE NEXT RUN
         G3S=G3S+DG3S
         FLAGP=1.0D0
C...
C...  IF (1 - G2(INF)) IS NEGATIVE, DECREASE G3(0)
      ELSE IF(FUNC(G2).LT.0.0D0)THEN
C...
C...     IF A SMALLER INCREMENT IN G3(0) IS REQUIRED, HALVE
C...     DG3S
         IF(FLAGP.GE.0.0D0)DG3S=DG3S/2.0D0
C...
C...     INCREMENT G3(0) NEGATIVELY FOR THE NEXT RUN
         G3S=G3S-DG3S
         FLAGM=1.0D0
      END IF
      RETURN
      END
      DOUBLE PRECISION FUNCTION FUNC(G2)
C...
C...  FUNCTION FUNC DEFINES THE FUNCTION FOR WHICH A ZERO IS TO BE
C...  COMPUTED
C...
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      FUNC=1.0D0-G2
      RETURN
      END
