--- a/rkmain.f90	Sun Jul 21 03:56:49 2024 +0900
+++ b/rkmain.f90	Sun Jul 21 04:11:29 2024 +0900
@@ -1,13 +1,12 @@
 program rkmain
   use rkmod
   implicit none
-  real(kind=dpkind),dimension(3) :: x
-  real(kind=dpkind) :: max_t
-  real(kind=dpkind),allocatable,dimension(:,:) :: r
+  real(kind=dpkind),dimension(3) :: x = [0.5, 0.5, 0.5]
+  integer :: max_i = 5000
+  real(kind=dpkind), allocatable, dimension(:, :) :: r
 
-  x = [0.5d0,0.5d0,0.5d0]
-  max_t = 5d1
-  call calc(x, max_t, r)
+  allocate(r(max_i, 4))
+  call calc(x, max_i, r)
   print *, r
 end program
 

--- a/rkmod.f90	Sun Jul 21 03:56:49 2024 +0900
+++ b/rkmod.f90	Sun Jul 21 04:11:29 2024 +0900
@@ -16,21 +16,21 @@
       rhside(3)=wk(1)*wk(2)-c*wk(3)
     end function
     !Rから呼び出すサブルーチン
-    subroutine calc(x, max_t, r)
+    subroutine calc(x, max_i, 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
+      integer, intent(in) ::  max_i
+      real(kind=dpkind) :: max_t
+      real(kind=dpkind),dimension(max_i, 4),intent(out) :: r
       real(kind=dpkind) :: zero = 0
-      integer i, max_i
+      integer i
       !微分方程式の右辺を計算する関数、解法、時間刻みの指定
       call rktype%setup(rhside,"rk4",1d-5) 
-      max_i = int(max_t / 0.01)
-      allocate(r(0:max_i, 4))
+      max_t = max_i * 0.01
       !時刻max_t秒まで、0.01秒毎に配列に記録。
       !ただし、内部的には刻み幅1e-5秒としている。
-      r(0, :) = [zero, x] ! 開始点
-      do i = 1, max_i
+      r(1, :) = [zero, x] ! 開始点
+      do i = 2, max_i
         call rktype%solver(rktype%time+1e-2,x) !積分実行
         r(i, :) = [rktype%time, x]
       end do