IMPLICIT NONE
      INTEGER,PARAMETER :: N=4,IMAX=10000000
      INTEGER :: I
      REAL(8) X(N),W(N*3),T,H,E
      REAL RX,RY
      EXTERNAL SUB

      CALL DCL

      T=0
      X(1)=3
      X(3)=0
      X(2)=0.3D0
      X(4)=0.2D0

      H=0.01D0
      
      DO I=1,IMAX
        CALL RK4(N,H,T,X,W,SUB)
        CALL ENE(X,E)
        IF(MOD(I,1000).EQ.0) THEN
          print *, I,E
          RX=X(1)
          RY=X(3)
          CALL SGPMZU(1,RX,RY,3,1,0.01)
        END IF
      END DO

      CALL GRCLS

      END
!-------------------------------------------------------------    
      SUBROUTINE DCL

      CALL SWPSET('IWIDTH',  400)
      CALL SWPSET('IHEIGHT', 400)
      CALL SWPSET('LWAIT',.FALSE.)
!      CALL SWPSET('LALT',.TRUE.)
      CALL SGPSET('LCORNER',.FALSE.)
      CALL GROPN(1)
      CALL UZPSET('INNER',-1)
      CALL SLRAT(1.0,1.0)
      CALL GLPSET('LMISS',.TRUE.)
      CALL SGPSET('LCNTL',.TRUE.)
      CALL SLRAT(1.0,1.0)
      
      CALL GRFRM
      CALL SGSVPT(0.2,0.9,0.2,0.9)
      CALL SGSWND(-1.0,4.0,-1.0,4.0)
      CALL SGSTRN(1)
      CALL SGSTRF
      CALL UXAXDV('B',0.5,1.0)
      CALL UXAXDV('T',0.5,1.0)
      CALL UYAXDV('L',0.5,1.0)
      CALL UYAXDV('R',0.5,1.0)
      CALL UXSTTL('B','X',0.0)
      CALL UYSTTL('L','Y',0.0)      
        
      END
!-------------------------------------------------------------    
      SUBROUTINE SUB(T,X,DX)

      IMPLICIT NONE
      INTEGER,PARAMETER :: N=4
      REAL(8) :: X(N),DX(N),T,R

      R=SQRT(X(1)*X(1)+X(3)*X(3))

      DX(1)=X(2)
      DX(3)=X(4)
      DX(2)=-X(1)/R**3
      DX(4)=-X(3)/R**3

      END
!-------------------------------------------------------------    
      SUBROUTINE RK4(N,H,T,X,W,SUB)

      IMPLICIT NONE
      REAL(8) :: X(N),W(N,3),T,H
      INTEGER :: I,N
      EXTERNAL SUB

      CALL SUB(T,X,W(1,1))
      DO I=1,N
        W(I,3)=X(I)+(H/2)*W(I,1)
      END DO

      T=T+H/2
      CALL SUB(T,W(1,3),W(1,2))
      DO I=1,N
        W(I,3)=X(I)+(H/2)*W(I,2)
        W(I,1)=W(I,1)+2*W(I,2)
      END DO

      CALL SUB(T,W(1,3),W(1,2))
      DO I=1,N
        W(I,3)=X(I)+H*W(I,2)
        W(I,1)=W(I,1)+2*W(I,2)
      END DO

      T=T+H/2
      CALL SUB(T,W(1,3),W(1,2))
      DO I=1,N
        X(I)=X(I)+(H/6)*(W(I,1)+W(I,2))
      END DO

      END
!-------------------------------------------------------------
      SUBROUTINE ENE(X,E)

      IMPLICIT REAL*8(A-H,O-Z)
      INTEGER,PARAMETER :: N=4
      REAL(8) :: X(N),R,E

      R=SQRT(X(1)*X(1)+X(3)*X(3))

      E=0.5D0*(X(2)*X(2)+X(4)*X(4))-1/R

      END