c  18.2$B85O"N)HyJ,J}Dx<0$N2rK!(B
c     $B?9@n(B $BLwBg(B , 22980138
c
      program ex18
c
      implicit none
      integer :: n 
      real :: x0,y0,z0,h
c
      n=10
      h=0.1e0
      x0=0.0e0
      y0=1.0e0
      z0=2.0e0
 
      write(6,*)'18. 2$B85O"N)HyJ,J}Dx<0$N2rK!(B'
      write(6,*)''
      call Runge_Kutta18  (x0,y0,z0,h,n)   

      stop 
      end
      
c-----$B%k%s%2!&%/%C%?K!$N%5%V%k!<%A%s(B
      subroutine Runge_Kutta18(x0,y0,z0,h,n)      
      implicit none
      integer :: n,i
      real :: h,xj,yj1,yj2,zj1,zj2,x0,y0,z0,f1,f2,
     -        k1,k2,k3,k4,m1,m2,m3,m4
      
      xj =x0
      yj1=y0
      zj1=z0
      write(6,200)
      write(6,100) 0.0,1.0e0,2.0e0
      do i=1,n
         k1=h*f1(xj,yj1,zj1)
         m1=h*f2(xj,yj1,zj1)
         k2=h*f1(xj+(h/2.0e0),yj1+(k1/2.0e0),zj1+(m1/2.0e0))
         m2=h*f2(xj+(h/2.0e0),yj1+(k1/2.0e0),zj1+(m1/2.0e0))
         k3=h*f1(xj+(h/2.0e0),yj1+(k2/2.0e0),zj1+(m2/2.0e0))
         m3=h*f2(xj+(h/2.0e0),yj1+(k2/2.0e0),zj1+(m2/2.0e0))
         k4=h*f1(xj+h,yj1+k3,zj1+m3)
         m4=h*f2(xj+h,yj1+k3,zj1+m3)

         yj2=yj1+(k1+2.0e0*k2+2.0e0*k3+k4)/6.0e0
         zj2=zj1+(m1+2.0e0*m2+2.0e0*m3+m4)/6.0e0
         xj =xj+h
         write(6,100) xj,yj2,zj2
         yj1=yj2
         zj1=zj2
      enddo
      
      return
  100 format(1x,f4.1,1x,e18.9,1x,e18.9) 
  200 format(4x,'x',4x,'y',18x,'z') 
      end

c-----$B4X?t(B f1=3y-2z $B$N%U%!%s%/%7%g%sJ8(B
      function f1(x,y,z)
      implicit none
      real :: f1,x,y,z

      f1=3.0e0*y-2.0e0*z

      return
      end

c-----$B4X?t(B f2=2y-z $B$N%U%!%s%/%7%g%sJ8(B
      function f2(x,y,z)
      implicit none
      real :: f2,x,y,z

      f2=2.0e0*y-z

      return
      end
