program rkmain
  use rungekutta !ルンゲクッタ解法をまとめたモジュール
  implicit none
  real(kind=dpkind),dimension(3) :: x !状態変数
  type(rk) :: rktype  !ルンゲクッタ法ソルバに関連する型
  real(kind=dpkind),allocatable,dimension(:,:) :: r
  real(kind=dpkind) :: max_t

  !初期値の指定
  x = [0.5e0,0.5e0,0.5e0]
  max_t = 5e1
  call calc(x, max_t, r)
  print *, r;

  !========================================
  contains
    !微分方程式の右辺を定義
    function rhside(self,time,wk)
      class(rk),intent(inout) :: self !parameter
      real(kind=dpkind) :: time
      real(kind=dpkind),dimension(:) :: wk
      real(kind=dpkind),dimension(size(wk)) :: rhside
      real(kind=dpkind),parameter :: a=10e0,b=28e0,c=8e0/3e0
      rhside(1)=a*(wk(2)-wk(1))
      rhside(2)=b*wk(1)-wk(2)-wk(1)*wk(3)
      rhside(3)=wk(1)*wk(2)-c*wk(3)
    end function
    !Rから呼び出すサブルーチン
    subroutine calc(x, max_t, r)
      implicit none
      real(kind=dpkind),dimension(3),intent(in) :: x
      real(kind=dpkind),intent(in) :: max_t
      real(kind=dpkind),allocatable,dimension(:,:),intent(out) :: r
      real(kind=dpkind) :: zero = 0
      integer i, max_i
      !微分方程式の右辺を計算する関数、解法、時間刻みの指定
      call rktype%setup(rhside,"rk4",1d-5) 
      max_i = int(max_t / 0.01)
      allocate(r(0:max_i, 4))
      !時刻max_t秒まで、0.01秒毎に配列に記録。
      !ただし、内部的には刻み幅1e-5秒としている。
      r(0, :) = [zero, x] ! 開始点
      do i = 1, max_i
        call rktype%solver(rktype%time+1e-2,x) !積分実行
        r(i, :) = [rktype%time, x]
      end do
    end subroutine
end program