!
! dU/dt = i omega U
! omega=pi
! U(0)=1.0
! 3段3位ルンゲ-クッタ公式
! ルンゲクッタ公式
!
! 2011/10/20? 荻原弘尭 作成
      program runge

      implicit none

      complex(8) :: kai1
      complex(8) :: kai2
      complex(8) :: kai3
      complex(8) :: kai4

      complex(8), allocatable :: U_kai(:)
      complex(8), parameter :: ei=(0.0d0, 1.0d0)
      real(8) :: omega
      real(8) :: pi
      real(8),allocatable :: time(:)
      real(8) :: dt
      real(8) :: t_max
      real(8) :: t_min
      real, allocatable :: U_r(:)
      real, allocatable :: time_r(:)
      integer :: i_l
      integer :: step
      integer :: IWS


      pi = 4.0d0*atan(1.0d0)
      omega=pi
      t_max=2.0d0*2.0d0*pi
      t_min=0.0d0

      write(*,*) "Enter in step"

      read(*,*) step

      allocate(U_kai(0:step), time(0:step),U_r(0:step),time_r(0:step))

      dt = (t_max-t_min)/dble(step)

      U_kai(0)=1.0d0
      time(0)=0.0d0
      do i_l = 1, step

      time(i_l) = t_min + dt*dble(i_l)

      kai1 = dt*omega*ei*U_kai(i_l-1)
      kai2 = dt*omega*ei*0.5d0*kai1
      kai3 = dt*omega*ei*0.5d0*kai2
      kai4 = dt*omega*ei*kai3

      U_kai(i_l) = &
        &  U_kai(i_l-1) + 1.0d0/6.0d0*( &
        & kai1 + 2.0d0*kai2+2.0d0*kai3 + kai4)

      write(*,*) "回数=",i_l,"U_kai=",real(U_kai(i_l))

      end do
      U_r(:)=real(U_kai(:))
      time_r(:)=real(time(:))

      call sgpwsn
      read(*,*) IWS
      call gropn(IWS)
      call grfrm
      call usspnt(step+1,time_r,U_r)
      call uspfit
      call grstrf
      call usdaxs
      call uulin(step+1,time_r,U_r)
      call grcls

      end program runge     
