      PROGRAM FIG4P2A
C...
C...  PROGRAM FIG4P2A CALCULATES AND PLOTS THE STABILITY DIAGRAM, I.E.,
C...  IM(LAMBDA*H) VS RE(LAMBDA*H) FOR THE BDF FORMULAS OF ORDERS ONE TO
C...  THREE, AS DESCRIBED IN GEAR, C. W., NUMERICAL INITIAL VALUE PROB-
C...  LEMS IN ORDINARY DIFFERENTIAL EQUATIONS, PRENTICE-HALL, ENGLEWOOD 
C...  CLIFFS, NJ, 1971, P217.  THUS, IT REPRODUCES FIGURE 11.6, PAGE 215,
C...  WHEN THE COEFFICIENTS OF TABLE 11.1 FOR ORDERS ONE TO THREE ARE
C...  USED, AS PLOTTED IN FIGURE 4.1A OF THE NUMERICAL METHOD OF LINES.
C...
C...  DIMENSION THE ARRAYS FOR THE BDF FORMULAS
      REAL ALPHA(6,6), BETA0(6)
C...
C...  OPEN A FILE FOR PLOTTING VIA TOP DRAWER
         OPEN(1,FILE='T.TOP',STATUS='NEW')
C...
C...  SCALE THE AXES OF THE PLOT
         WRITE(1,6)
6        FORMAT(' SET LIMITS X FROM -2 TO 6 Y FROM -4 TO 4',/,
     1          ' SET FONT DUPLEX')
         WRITE(1,7)
7        FORMAT(' SET SIZE 6 BY 6')
C...
C...  SET THE COEFFICIENTS OF THE BDF FORMULAS (FROM TABLE 11.1, P217,
C...  GEAR)
C...
C...     FIRST ORDER METHOD
         ALPHA(1,1)=1.0
           BETA0(1)=1.0
C...
C...     SECOND ORDER METHOD
         ALPHA(1,2)= 4.0/3.0
         ALPHA(2,2)=-1.0/3.0
           BETA0(2)= 2.0/3.0
C...
C...     THIRD ORDER METHOD
         ALPHA(1,3)= 18.0/11.0
         ALPHA(2,3)= -9.0/11.0
         ALPHA(3,3)=  2.0/11.0
           BETA0(3)=  6.0/11.0
C...
C...  STEP THROUGH ORDERS 1 TO 3
         DO 10 NORD=1,3
C...
C...  INCREMENT IN THETA IS 2*PI/300
      DTHETA=2.0*4.0*ATAN(1.0)/300.0
C...
C...  STEP THROUGH THE VALUES OF THETA FROM O TO 2*PI
      DO 4 I=1,301
      THETA=FLOAT(I-1)*DTHETA
C...
C...  COMPUTE THE REAL AND IMAGINARY PARTS OF H*LAMBDA
         HLREAL=1.0
         HLIMAG=0.0
         DO 2 IT=1,NORD
         HLREAL=HLREAL-ALPHA(IT,NORD)*COS(FLOAT(IT)*THETA)
         HLIMAG=HLIMAG+ALPHA(IT,NORD)*SIN(FLOAT(IT)*THETA)
C...
C...  PRINT INTERMEDIATE VALUES OF REAL AND IMAGINARY PARTS (TO FOLLOW
C...  THE CALCULATIONS IN DETAIL)
C        WRITE(*,1)IT,HLREAL,HLIMAG
1        FORMAT(' IT = ',I3,'  REAL = ',F10.4,'  IMAG = ',F10.4)
2        CONTINUE
C...
C...  DIVISION BY BETA0 GIVES THE FINAL REAL AND IMAGINARY PARTS
         HLREAL=HLREAL/BETA0(NORD)
         HLIMAG=HLIMAG/BETA0(NORD)
C...
C...  PRINT THE REAL AND IMAGINARY PARTS (OF H*LAMBDA)
         WRITE(*,3)THETA,HLREAL,HLIMAG
3        FORMAT('  THETA = ',F7.3,
     1          '   REAL = ',F7.3,
     2          '   IMAG = ',F7.3)
C...
C...  WRITE THE REAL AND IMAGINARY PARTS TO A FILE FOR TOP DRAWER 
C...  PLOTTING
         WRITE(1,5)HLREAL,HLIMAG
5        FORMAT(2F10.4)
C...
C...  GO ON TO THE NEXT VALUE OF THETA
4     CONTINUE
C...
C...  CONNECT THE POINTS
      WRITE(1,8)
8     FORMAT(' JOIN')
C...
C...  GO ON TO THE NEXT STABILITY BOUNDARY
10    CONTINUE
C...
C...  DRAW A VERTICAL LINE THROUGH (0,0)
         X=0.
         DO 12 I=1,11
         Y=-4.0+8.0*FLOAT(I-1)/10.0
         WRITE(1,5)X,Y
12       CONTINUE
         WRITE(1,11)
11       FORMAT(' JOIN')
C...
C...  DRAW A HORIZONTAL LINE THROUGH (0,0)
         Y=0.
         DO 13 I=1,11
         X=-2.0+8.0*FLOAT(I-1)/10.0
         WRITE(1,5)X,Y
13       CONTINUE
         WRITE(1,11)
      STOP
      END
