c  17.1$B3,>oHyJ,J}Dx<0$N=i4|CMLdBj$N(B4$B<oN`$N2rK!(B
c     $B?9@n(B $BLwBg(B , 22980138
c
      program ex17
c
      implicit none
      integer :: n 
      real :: x0,y0,h
c
      n=10
      h=0.1e0
      x0=1.0e0
      y0=1.0e0

      write(6,*)'17.1$B3,>oHyJ,J}Dx<0$N(B4$B<oN`$N2rK!(B'
      write(6,*)''
      write(6,*)'(1)$B%*%$%i!<K!(B'
      call Euler17(x0,y0,h,n)

      write(6,*)''
      write(6,*)'(2)$B%[%$%sK!(B'
      call Heun17(x0,y0,h,n)

      write(6,*)''
      write(6,*)'(3)$B2~NI%*%$%i!<K!(B'
      call Improve_Euler17(x0,y0,h,n)

      write(6,*)''
      write(6,*)'(4)$B%k%s%2!&%/%C%?K!(B'
      call Runge_Kutta17  (x0,y0,h,n)   

      stop 
      end
      
c-----$B%*%$%i!<K!$N%5%V%k!<%A%s(B
      subroutine Euler17(x0,y0,h,n)
      implicit none
      integer :: n,i
      real :: h,xj,yj1,yj2,x0,y0,f,e
      
      xj =x0
      yj1=y0
      write(6,200)
      e=abs(sqrt(1.0e0)-1.0e0)
      write(6,100) 1.0,1.0e0,sqrt(1.0e0),e
      do i=1,n
         yj2=yj1+h*f(xj,yj1)
         xj =xj+h
         e=abs(sqrt(xj)-yj2)
         write(6,100) xj,yj2,sqrt(xj),e
         yj1=yj2
      enddo

      return
  100 format(1x,f4.1,1x,e18.9,1x,e18.9,1x,e18.9) 
  200 format(4x,'x',4x,'y',18x,'sqrt(x)',12x,'$B8m:9(B') 
      end

c-----$B%[%$%sK!$N%5%V%k!<%A%s(B
      subroutine Heun17(x0,y0,h,n)      
      implicit none
      integer :: n,i
      real :: h,xj,yj1,yj2,x0,y0,f,k1,k2,e
      
      xj =x0
      yj1=y0
      write(6,200)
      e=abs(sqrt(1.0e0)-1.0e0)
      write(6,100) 1.0,1.0e0,sqrt(1.0e0),e
      do i=1,n
         k1=h*f(xj,yj1)
         k2=h*f(xj+h,yj1+h)
         yj2=yj1+(k1+k2)/2.0e0
         xj =xj+h
         e=abs(sqrt(xj)-yj2)
         write(6,100) xj,yj2,sqrt(xj),e
         yj1=yj2
      enddo
      
      return
  100 format(1x,f4.1,1x,e18.9,1x,e18.9,1x,e18.9) 
  200 format(4x,'x',4x,'y',18x,'sqrt(x)',12x,'$B8m:9(B') 
      end

c-----$B2~NI%*%$%i!<K!$N%5%V%k!<%A%s(B
      subroutine Improve_Euler17(x0,y0,h,n)      
      implicit none
      integer :: n,i
      real :: h,xj,yj1,yj2,x0,y0,f,k1,k2,e
      
      xj =x0
      yj1=y0
      write(6,200)
      e=abs(sqrt(1.0e0)-1.0e0)
      write(6,100) 1.0,1.0e0,sqrt(1.0e0),e
      do i=1,n
         k1=h*f(xj,yj1)
         k2=h*f(xj+(h/2.0e0),yj1+(k1/2.0e0))
         yj2=yj1+k2
         xj =xj+h
         e=abs(sqrt(xj)-yj2)
         write(6,100) xj,yj2,sqrt(xj),e
         yj1=yj2
      enddo
      
      return
  100 format(1x,f4.1,1x,e18.9,1x,e18.9,1x,e18.9) 
  200 format(4x,'x',4x,'y',18x,'sqrt(x)',12x,'$B8m:9(B') 
      end

c-----$B%k%s%2!&%/%C%?K!$N%5%V%k!<%A%s(B
      subroutine Runge_Kutta17(x0,y0,h,n)      
      implicit none
      integer :: n,i
      real :: h,xj,yj1,yj2,x0,y0,f,k1,k2,k3,k4,e
      
      xj =x0
      yj1=y0
      write(6,200)
      e=abs(sqrt(1.0e0)-1.0e0)
      write(6,100) 1.0,1.0e0,sqrt(1.0e0),e
      do i=1,n
         k1=h*f(xj,yj1)
         k2=h*f(xj+(h/2.0e0),yj1+(k1/2.0e0))
         k3=h*f(xj+(h/2.0e0),yj1+(k2/2.0e0))
         k4=h*f(xj+h,yj1+k3)
         yj2=yj1+(k1+2.0e0*k2+2.0e0*k3+k4)/6.0e0
         xj =xj+h
         e=abs(sqrt(xj)-yj2)
         write(6,100) xj,yj2,sqrt(xj),e
         yj1=yj2
      enddo
      
      return
  100 format(1x,f4.1,1x,e18.9,1x,e18.9,1x,e18.9) 
  200 format(4x,'x',4x,'y',18x,'sqrt(x)',12x,'$B8m:9(B') 
      end

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

      f=y/(2.0e0*x)

      return
      end
