      PROGRAM FIG4P1
C...
C...  PROGRAM FIG4P1 CALCULATES AND PLOTS THE STABILITY DIAGRAM, I.E.,
C...  IM(LAMBDA*H) VS RE(LAMBDA*H) FOR SOME EXPLICIT ONE STEP INTEGRA-
C...  TION METHODS, E.G., EXPLICIT EULER, CLASSICAL FOURTH ORDER RUNGE
C...  KUTTA, USING THE METHOD DESCRIBED BY GEAR, C. W., NUMERICAL INITIAL
C...  VALUE PROBLEMS IN ORDINARY DIFFERENTIAL EQUATIONS, PRENTICE-HALL,
C...  ENGLEWOOD CLIFF, NJ, 1971, P41.  THUS, IT REPRODUCES FIGURE 2.1,
C...  PAGE 41, WHEN THE CLASSICAL FOURTH ORDER RUNGE KUTTA METHOD IS
C...  SELECTED, AS PLOTTED IN FIGURE 4.1 OF THE NUMERICAL METHOD OF LINES.
C...
C...  ORDER OF SYSTEM
      PARAMETER (NDEG=4)
C...
C...  ARRAYS FOR COEFFICIENTS AND ZEROS (ROOTS)
      COMPLEX COEFF(NDEG+1), ZERO(NDEG)
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 -3 TO 1 Y FROM -3.5 TO 3.5',/,
     1          ' SET FONT DUPLEX')
         WRITE(1,7)
7        FORMAT(' SET SIZE 6 BY 6')
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 2 I=1,301
      THETA=FLOAT(I-1)*DTHETA
C...
C...  SET COEFFICIENTS FOR THE CLASSICAL FOURTH RUNGE KUTTA METHOD (FOR 
C...  THE EXPLICIT EULER METHOD, SET ONLY COEFF(1) AND COEFF(2), AND SET
C...  NDEG = 1 AT THE BEGINNING)
      COEFF(1)=CMPLX(1.0-COS(THETA),-SIN(THETA))
      COEFF(2)=CMPLX(1.0,0.0)
      COEFF(3)=CMPLX(1.0/2.0,0.)
      COEFF(4)=CMPLX(1.0/6.0,0.)
      COEFF(5)=CMPLX(1.0/24.0,0.)
C...
C...  CALCULATE ROOTS
      CALL ZPOCC(NDEG,COEFF,ZERO)
C...
C...  PRINT INDEX FOR THETA ON THE SCREEN (TO OBSERVE THE PROGRESS OF 
C...  THE CALCULATION) 
      WRITE(*,10)I
10    FORMAT(I4)
C...
C...  PRINT ROOTS ON THE SCREEN (TO EXAMINE THEIR VALUES IN DETAIL)
      DO 8 IR=1,NDEG
         AR= REAL(ZERO(IR))
         AI=AIMAG(ZERO(IR))
C        WRITE(*,1)I,AR,AI
1        FORMAT(' I = ',I3,'  REAL = ',F10.4,'  IMAG = ',F10.4)
8     CONTINUE
C...
C...  WRITE THE ROOTS TO A FILE FOR TOP DRAWER PLOTTING
      DO 9 IR=1,NDEG
         AR= REAL(ZERO(IR))
         AI=AIMAG(ZERO(IR))
         WRITE(1,3)AR,AI
3        FORMAT(2F10.4)
9     CONTINUE
C...
C...  GO ON TO THE NEXT VALUE OF THETA
2     CONTINUE
C...
C...  CONNECT THE POINTS
      WRITE(1,5)
5     FORMAT(' PLOT')
C...
C...  DRAW A VERTICAL LINE THROUGH (0,0)
         X=0.
         DO 12 I=1,11
         Y=-3.5+7.0*FLOAT(I-1)/10.0
         WRITE(1,3)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=-3.0+4.0*FLOAT(I-1)/10.0
         WRITE(1,3)X,Y
13       CONTINUE
         WRITE(1,11)
        STOP
      END
