C R.C.
      PROGRAM OSCIATION
      PARAMETER( NMAX=10000 )
      REAL Y(0:NMAX), U(0:NMAX), X(0:NMAX)
      REAL K1(0:NMAX), K2(0:NMAX), K3(0:NMAX), K4(0:NMAX)
      REAL J1(0:NMAX), J2(0:NMAX), J3(0:NMAX), J4(0:NMAX)
      REAL A, B, gamma, kappa
C 配列を宣言するとき、配列の長さは変数で指定できない。ここでの変数は NMAX.
C そのため、配列の長さを変数で指定したい場合は、PARAMETER宣言を行う.
      REAL DX
      PI=3.14159265
*-- 減衰係数および振動数を指定する ---
      WRITE(*,*) "INPUT THE DAMPED AND FREQUENCY."
      READ(*,*) kappa, B
      A=B**2
      gamma=2.0*kappa
      XMAX=10.0
      XMIN=0.0
      R=XMAX-XMIN
      DX=R/NMAX
C 以下 4 行が境界条件
      Y(0)=1.0
      U(0)=10.0
      X(0)=0.0
      DO 10 K=0,NMAX-1
        K1(K)=DX*(U(K))
        J1(K)=DX*(A*(-SIN(Y(K)))-gamma*U(K))
        K2(K)=DX*(U(K)+0.5*J1(K))
        J2(K)=DX*(A*(-SIN(Y(K)+0.5*K1(K)))-gamma*(U(K)+0.5*J1(K)))
        K3(K)=DX*(U(K)+0.5*J2(K))
        J3(K)=DX*(A*(-SIN(Y(K)+0.5*K2(K)))-gamma*(U(K)+0.5*J2(K)))
        K4(K)=DX*(U(K)+J3(K))
        J3(K)=DX*(A*(-SIN(Y(K)+K3(K)))-gamma*(U(K)+J3(K)))
        X(K+1)=X(K)+DX
        Y(K+1)=Y(K)+(1.0/6.0)*(K1(K)+2*K2(K)+2*K3(K)+K4(K))
        U(K+1)=U(K)+(1.0/6.0)*(J1(K)+2*J2(K)+2*J3(K)+J4(K))
  10  CONTINUE
*-- グラフの呼び出し ---
      WRITE(*,*) ' WORKSTATION ID (I)  ? ;'
      CALL SGPWSN
      READ (*,*) IWS
      CALL GROPN( IWS )
      CALL GRFRM
      CALL USSPNT( NMAX+1, X, Y )
      CALL GRSWND( XMIN, XMAX, 0.0, 10.0 )
      CALL GRSVPT(  0.2,  0.8,  0.2,  0.8 )
      CALL USPFIT
      CALL GRSTRF
      CALL USSTTL( '', '', '', '' )
      CALL USDAXS
      CALL UULIN( NMAX+1, X, Y )
      CALL GRCLS
      STOP
      END