--- a/rkmain.f90	Tue Jul 23 02:30:57 2024 +0900
+++ b/rkmain.f90	Sun Jul 21 02:08:58 2024 +0900
@@ -3,18 +3,15 @@
   implicit none
   real(kind=dpkind),dimension(3) :: x !状態変数
   type(rk) :: rktype  !ルンゲクッタ法ソルバに関連する型
+  real(kind=dpkind),allocatable,dimension(:,:) :: r
+  real(kind=dpkind) :: max_t
 
-  !微分方程式の右辺を計算する関数、解法、時間刻みの指定
-  call rktype%setup(rhside,"rk4",1d-5) 
   !初期値の指定
-  x=[0.5e0,0.5e0,0.5e0]
-  !時刻0?50秒まで、0.01秒毎に標準出力へ書き出し。
-  !ただし、内部的には刻み幅1e-5秒としている。
-  do 
-    if(rktype%time>5e1)exit
-    call rktype%solver(rktype%time+1e-2,x) !積分実行
-    print*,rktype%time,x
-  enddo
+  x = [0.5e0,0.5e0,0.5e0]
+  max_t = 5e1
+  call calc(x, max_t, r)
+  print *, r;
+
   !========================================
   contains
     !微分方程式の右辺を定義
@@ -28,4 +25,25 @@
       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
+