program ex042
C
C       This program solves the differential equation
C           dy/dx = y/2x
C       by using the modified Euler method.
C
      implicit real*8(a-h,o-z)
C
      xi = 1.0d0            ! Initial value of x
      xf = 2.0d0            ! Final   value of x
      yi = 1.0d0            ! Initial value of y
C
      open(16,file='ex042.dat')     ! open(unit,file='filename')
C
      do n = 10, 100, 10
          dx = (xf-xi)/n    ! Mesh size
          x = xi            ! Initial Conditions
          y = yi            ! Initial Conditions
          do k = 1, n
              dy1 = f(x,y)*dx
              x1 = x + dx /2    ! x' = x(k-1) + dx / 2
              y1 = y + dy1/2    ! y' = y(k-1) + dx * f(x(k-1), y(k-1)) / 2
c
              dy = f(x1,y1)*dx
              x = x + dx        ! x(k) = x(k-1) + dx
              y = y + dy        ! y(k) = y(k-1) + dx * f(x(k-1), y(k-1))
          enddo
          write( *,*) dx,abs(y-sqrt(x))
          write(16,*) dx,abs(y-sqrt(x))
      end do
C
      stop
      end
c **********************************************************************
      function f(x,y)
c **********************************************************************
      implicit real*8 (a-h,o-z)
      f = 0.5d0*y/x
      end