      PROGRAM EULER1
*
*     Double precision coding is used
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
*
*     Open a file for output
      OPEN(1,FILE='pend1.out')
*
*     Problem parameters
*
*        Pendulum length (m)
         XL=30.0D0
*
*        Acceleration of gravity (m/s**2)
         G=9.8D0
*
*     Step through a series of Euler integrations
*     (with varying step size, h)
      DO IH=1,4
         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
*
*        Begin integration (set initial conditions)
*
*        Initial time (sec)
         T=0.0D0
*
*        Initial angle (radians)
*           Linear problem
            Y3=1.0D0
            Y30=Y3
*
*        Initial velocity (radians/sec)
*           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(//,7X,'h',7X,'t',8X,'y3',7X,'y3e',6X,'error')
            ERROR=0.0D0
            WRITE(*,1)H,T,Y3,Y30,ERROR
            WRITE(1,1)H,T,Y3,Y30,ERROR
1           FORMAT(F8.4,F8.2,2F10.4,F11.6)
*
*           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
*
*           Exact solution and error
            Y3E=Y30*DCOS(DSQRT(G/XL)*T)
            ERROR=Y3-Y3E
*
*           Write numerical and exact solutions, error
            WRITE(*,1)H,T,Y3,Y3E,ERROR
            WRITE(1,1)H,T,Y3,Y3E,ERROR
            END IF
*
*        Derivative at base point
*           Linear problem
            DY3DT=Y4
            DY4DT=-(G/XL)*Y3
*
*        Euler step
         T=DFLOAT(I+1)*H
*
*           Linear problem
            Y3=Y3+DY3DT*H
            Y4=Y4+DY4DT*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,
     +  'Euler''s method is first order correct,',/,3X,
     +  'i.e., O(h)')
*
*     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
