program main
      implicit none

      interface
        function dy(t,y)
          real,intent(in) ::t,y
          real :: dy
        end function dy
      end interface

      real :: t0=0.0,y0=3.0
      integer :: n=200
      real :: tr,yr     !runge-kutta
      real :: te,ye     !euler
      integer :: i      !do-loop
      real :: k1,k2,k3,k4 !rk-method use
      real :: h
        h=(3.0-1.0)/real(n)/2*3


      open(unit=15,file='final',status='unknown')
        write(15,"('#',6x,'t',14x,'y')")
        write(*,"('#',6x,'t',14x,'y')")


!euler-----------------------------------------------------------

        te = 0.0
        ye = 0.0
        do i=1,n+1
          ye = y0 + h*dy(te,ye)
          te = t0 + (i-1.0)*h
!           if(te >= 1.0) write(*,*) te,ye
          y0 = ye
        end do

!runge-kutta----------------------------------------------------
        tr = 0.0
        yr = 0.0
        k1 = 0.0
        k2 = 0.0
        k3 = 0.0
        k4 = 0.0
        do i=1,n+1
          tr = t0 + (i-1) * h
          yr = y0 + (h/6.0) * ( k1 + 2*k2 + 2*k3 + k4 )
          k1 = dy(tr,y0)
          k2 = dy(tr+h/2.0 , y0+h/2.0*k1)
          k3 = dy(tr+h/2.0 , y0+h/2.0*k2)
          k4 = dy(tr+h , y0+h*k3)
        do i=1,n+1
          tr = t0 + (i-1) * h
          yr = y0 + (h/6.0) * ( k1 + 2*k2 + 2*k3 + k4 )
          k1 = dy(tr,y0)
          k2 = dy(tr+h/2.0 , y0+h/2.0*k1)
          k3 = dy(tr+h/2.0 , y0+h/2.0*k2)
          k4 = dy(tr+h , y0+h*k3)
            if (tr >= 1.0) write(*,*) tr,yr
          y0 = yr
        end do
      close(unit=15)

      end program main
!================================================================
      real function dy(t,y)
      implicit none
      real :: t,y
        dy = 2.0*t*y + t
      return
      end function dy