program list8_3
  !ordinary differential equation main routine
  !runge-kutta method
  implicit none
  integer ni,i
  integer:: ind=0
  real(8) xu,xl,d,S(20),x
  data xu,xl,d,S(1),ni/2.0,1.0,0.5,1.0,1/
  write(*,*)' *** ordinary differential equation ***'
  write(*,*)'        (runge-kutta method)'
  write(*,*)
  write(*,200) xu,xl,d,S(1),ni
  write(*,*)
  call rugkut(xu,xl,d,S,ni,ind)
  if(ind.ne.0)then
     write(*,*)'     ind=',ind
     stop
  endif
  write(*,*)' * solution *'    
  i=1
  x=xl
10 write(*,220) x,S(i)
  x=x+d
  i=i+1
  if(x.le.xu)goto 10
  stop
200 format(1H ,'upper limit=',f5.1/'lower limit=',f5.1/'sampling interval=',f6.2/'initial condition=',f5.1/'subdivision number=',i3)

220 format(1H ,'x=',f5.1,3x,'y=',f8.5)
end program list8_3


!function f
real(8) function F(x,y)
  implicit none
  real(8):: F
  real(8) x,y,z
  !external F
  z = - x * y + y**2  !定義式
  F = z

  return
end function F


!ordinary differential equation subroutine
!runge-kutta method
subroutine rugkut(xu,xl,d,S,ni,ind)
  implicit none
  real(8) S(20),x,y,z,w,h,s1,s2,s3,s4,xu,xl,d
  real(8)::F
  integer l,k,ni,n
  integer ind
  if(n.gt.20)then
     ind=54
     return
  endif
  x=xl
  y=S(1)
  z=xl
  h=d/real(ni)
  l=2
10 if(l.gt.20)then
     ind=55
     return
  endif
  do k=1,ni
     w=y
     s1=h*F(x,w)
     w=y+0.5*s1
     s2=h*F(x+0.5*h,w)
     w=y+0.5*s2
     s3=h*F(x+0.5*h,w)
     w=y+s3
     s4=h*F(x+h,w)
     y=y+(s1+2.0*s2+2.0*s3+s4)/6.0
     x=x+h
  enddo
     S(l)=y
     z=z+d
     if(z.ge.xu)return
     x=z
     l=l+1
     goto 10
   end subroutine rugkut



  *** ordinary differential equation ***
         (runge-kutta method)

upper limit=  2.0
lower limit=  1.0
sampling interval=  0.50
initial condition=  1.0
subdivision number=  1

  * solution *
 x=  1.0   y= 1.00000
 x=  1.5   y= 0.86352
 x=  2.0   y= 0.51018