--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/rkmain.f90	Tue Jul 23 02:30:57 2024 +0900
@@ -0,0 +1,31 @@
+program rkmain
+  use rungekutta !ルンゲクッタ解法をまとめたモジュール
+  implicit none
+  real(kind=dpkind),dimension(3) :: x !状態変数
+  type(rk) :: rktype  !ルンゲクッタ法ソルバに関連する型
+
+  !微分方程式の右辺を計算する関数、解法、時間刻みの指定
+  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
+  !========================================
+  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
+end program