PROGRAM EULER5
*
*     Double precision coding is used
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
*
*     Open a file for output
      OPEN(1,FILE='pend5.out')
*
*     Problem parameters
*
*        Pendulum length (m)
         XL=30.0D0
*
*        Acceleration of gravity (m/s**2)
         G=9.8D0
*
*     Step through a series of modified Euler integrations
*     (with varying step size, h)
      DO IH=1,5
         IF(IH.EQ.1)H=1.0D0
         IF(IH.EQ.2)H=0.1D0
         IF(IH.EQ.3)H=0.01D0
         IF(IH.EQ.4)H=0.001D0
         IF(IH.EQ.5)H=0.0001D0
*
*        Begin integration (set initial conditions)
*
*        Initial time (sec)
         T=0.0D0
*
*        Initial angle (radians)

*           Nonlinear problem
            Y1=1.0D0
*
*           Linear problem
            Y3=1.0D0
*
*        Initial velocity (radians/sec)
*
*           Nonlinear problem
            Y2=0.0D0
*
*           Linear problem
            Y4=0.0D0
*
*        Final t (sec)
         TF=5.0D0
*
*        Number of steps for 0 le t le TF
         NSTEPS=INT(TF/H)
*
*        Perform integration
         DO I=0,NSTEPS
*
*           Write numerical solutions
*
*           t=0
            IF(T.LT.0.00001D0)THEN
            WRITE(*,2)
            WRITE(1,2)
2           FORMAT(//,9X,'h',7X,'t',13X,'y1',13X,'y3')
            WRITE(*,1)H,T,Y1,Y3
            WRITE(1,1)H,T,Y1,Y3
1           FORMAT(F10.4,F8.2,2F15.9)
*
*           t = 1, 2, 3, 4, 5
            ELSE IF (
     +      ((T.GT.0.99999D0).AND.(T.LT.1.00001D0)).OR.
     +      ((T.GT.1.99999D0).AND.(T.LT.2.00001D0)).OR.
     +      ((T.GT.2.99999D0).AND.(T.LT.3.00001D0)).OR.
     +      ((T.GT.3.99999D0).AND.(T.LT.4.00001D0)).OR.
     +      ((T.GT.4.99999D0).AND.(T.LT.5.00001D0)))THEN
            WRITE(*,1)H,T,Y1,Y3
            WRITE(1,1)H,T,Y1,Y3
            END IF
*
*        Derivative at base point
*
*           Nonlinear problem
            Y1B=Y1
            Y2B=Y2
            DY1DTB=Y2B
            DY2DTB=-(G/XL)*DSIN(Y1B)
*
*           Linear problem
            Y3B=Y3
            Y4B=Y4
            DY3DTB=Y4B
            DY4DTB=-(G/XL)*Y3B
*
*        Euler step
         T=DFLOAT(I+1)*H
*
*           Nonlinear problem
            Y1=Y1B+DY1DTB*H
            Y2=Y2B+DY2DTB*H
*
*           Linear problem
            Y3=Y3B+DY3DTB*H
            Y4=Y4B+DY4DTB*H
*
*        Derivative at advanced points
*
*           Nonlinear problem
            DY1DT=Y2
            DY2DT=-(G/XL)*DSIN(Y1)
*
*           Linear problem
            DY3DT=Y4
            DY4DT=-(G/XL)*Y3
*
*        Modified Euler step
*
*           Nonlinear problem
            Y1=Y1B+(DY1DTB+DY1DT)/2.0D0*H
            Y2=Y2B+(DY2DTB+DY2DT)/2.0D0*H
*
*           Linear problem
            Y3=Y3B+(DY3DTB+DY3DT)/2.0D0*H
            Y4=Y4B+(DY4DTB+DY4DT)/2.0D0*H
*
*        Continue Euler integration
         END DO
*
*     Integration complete; next integration step size
      WRITE(*,3)NSTEPS
      WRITE(1,3)NSTEPS
3     FORMAT(/,3X,'NSTEPS = ',I6)
      PAUSE
      END DO
*
*     Calculation for all integration steps completed
      WRITE(*,4)
      WRITE(1,4)
4     FORMAT(///,'Conclusion:                               ',/,3X,
     +  'The modified Euler method can easily be programmed ',/,3X,
     +  'for systems of nonlinear ODEs for which analytical ',/,3X,
     +  'solutions and associated theorems are unavailable; ',/,3X,
     +  'thus, for realistic engineering problems, we must  ',/,3X,
     +  'use numerical methods.                            ',//,3X,
     +  'Plotting the solutions for the linear and nonlinear',/,3X,
     +  'problems would demonstrate the differences in the  ',/,3X,
     +  'solutions, e.g., the nonlinear problem would not   ',/,3X,
     +  'have pure linear harmonics, i.e., sine and cosine  ',/,3X,
     +  'solutions.                                         ')
*
*     End of calculations
      STOP
      END
      DOUBLE PRECISION FUNCTION DFLOAT(I)
*
*     Function DFLOAT converts a single precision integer to a
*     double precision float
*
      DFLOAT=DBLE(REAL(I))
      RETURN
      END