program main
      integer :: i,n
      real :: x,y,z,d,tmax
      real,dimension(:),allocatable :: px,py,pz
      intrinsic nint                        !実数を四捨五入する. 配列の大きさを求めるのに使用
      external Runge
      x=1.0; y=1.0; z=1.0; d=0.01; tmax=100.0  !(x,y,z)は初期条件, dは刻み幅, xmaxまで計算する
      n=nint(tmax/d)+1                      ! +1で初期条件の分も含めた配列の大きさに
      allocate(px(n),py(n),pz(n))           !配列の大きさを決定
      call runge(x,y,z,d,n,px,py,pz)      !関数の処理が終わると配列3つに結果が格納される
      open(9,file='rolents.data')
      do i=1,n
        write(9,*) px(i),py(i),pz(i)        ! t0+dからの位置(x,y,z)を書き出し
      end do
      endfile(9)
      deallocate(px,py,pz)
end program main

!---------関数を定義--------------------------
function f1(t,x,y,z) result(res)
      real :: t,x,y,z,res
      res=10.*y-10.*x
end function f1
!----------------------------------------------
function f2(t,x,y,z) result(res)
        real :: t,x,y,z,res
        res=28.*x-y-x*z
end function f2
!-------------------------------------------
function f3(t,x,y,z) result(res)
        real :: t,x,y,z,res
        res=-8./3.*z+x*y
end function f3

!---------サブルーチンを定義----------------
subroutine runge(x,y,z,d,n,px,py,pz)
      real,intent(in) :: x,y,z,d
      integer,intent(in) :: n     
      real,dimension(n),intent(out) :: px,py,pz
      integer :: i
      real :: t0,x0,y0,z0
      real :: d11,d12,d13,d14
      real :: d21,d22,d23,d24
      real :: d31,d32,d33,d34
      real,external :: f1,f2,f3
      x0=x; y0=y; z0=z; t0=0.0
      px(1)=x0; py(1)=y0; pz(1)=z0
      do i=2,n
        d11=d*f1(t0,x0,y0,z0)
        d21=d*f2(t0,x0,y0,z0)
        d31=d*f3(t0,x0,y0,z0)

        d12=d*f1(t0+d/2.,x0+d11/2.,y0+d21/2.,z0+d31/2.)
        d22=d*f2(t0+d/2.,x0+d11/2.,y0+d21/2.,z0+d31/2.)
        d32=d*f3(t0+d/2.,x0+d11/2.,y0+d21/2.,z0+d31/2.)

        d13=d*f1(t0+d/2.,x0+d12/2.,y0+d22/2.,z0+d32/2.)
        d23=d*f2(t0+d/2.,x0+d12/2.,y0+d22/2.,z0+d32/2.)
        d33=d*f3(t0+d/2.,x0+d12/2.,y0+d22/2.,z0+d32/2.)

        d14=d*f1(t0+d,x0+d13,y0+d23,z0+d33)
        d24=d*f2(t0+d,x0+d13,y0+d23,z0+d33)
        d34=d*f3(t0+d,x0+d13,y0+d23,z0+d33)

        t0=t0+d
        x0=x0+(d11+2.*d12+2.*d13+d14)/6.
        y0=y0+(d21+2.*d22+2.*d23+d24)/6.
        z0=z0+(d31+2.*d32+2.*d33+d34)/6.
        px(i)=x0
        py(i)=y0
        pz(i)=z0
      end do
end subroutine runge
