*-----------------------------------------------------------------------
*     4th order Runge-Kutta algorithm routine.
*                                                 Oct. 5, 1990  S.Sakai
*-----------------------------------------------------------------------
*     Copyright (C) 2000-2016 GFD Dennou Club. All rights reserved.
*-----------------------------------------------------------------------
      SUBROUTINE ODRK4(N, FCN, T, DT, X, DX, XOUT, WORK)
      DIMENSION  X(N), XOUT(N), DX(N), WORK(N,3)

      DTT = DT/2.
      TT  = T+DTT
      DO 10 I=1,N
        WORK(I,1) = X(I) + DTT*DX(I)
   10 CONTINUE
      CALL FCN(N, TT, WORK(1,1), WORK(1,2))

      DO 20 I=1,N
        WORK(I,1) = X(I) + DTT*WORK(I,2)
   20 CONTINUE
      CALL FCN(N, TT, WORK(1,1), WORK(1,3))

      TT = T+DT
      DO 30 I=1,N
        WORK(I,1) = X(I) + DT*WORK(I,3)
        WORK(I,3) = WORK(I,2) + WORK(I,3)
   30 CONTINUE
      CALL FCN(N, TT, WORK(1,1), WORK(1,2))

      DTT=DT/6.
      DO 40 I=1,N
        XOUT(I) = X(I) + DTT*(DX(I) + 2.*WORK(I,3) + WORK(I,2))
  40  CONTINUE

      RETURN
      END
