!**********************************************************************
 Program ode
!**********************************************************************
	implicit none
	integer::i,N
	real(8)::x,x0,y,y0,xmax
	real(8),allocatable::xval(:),yval(:)
	character(80)::FM,FM1
!----------------------------------------------------------------------

	interface 
	real(8) function f(x,y)
	real(8),intent(in)::x,y
	end function f
	end interface
!----------------------------------------------------------------------
!Kezdeti feltétel
    x0=0.d0 ; y0=1.d0
!Osztáspontok száma 
    N=5  
!A tartomány [x0,xmax] x0<xmax
    xmax=1.d0
! tömbök allokálása
	allocate(xval(N),yval(N))
! Call Runge-Kutta
	call RK4(x0,xmax,y0,N,xval,yval)
!  Print result
	FM='(15x,"x",20x,"y",/,5x,40("-"))'
	FM1='(5x,f15.10,5x,f15.10)'
  print FM
  Print FM1,(xval(i),yval(i),i=1,N)
  deallocate(xval,yval)
 end program ode
!**********************************************************************
 real(8) function f(x,y)
!**********************************************************************
	real(8),intent(in)::x,y
!----------------------------------------------------------------------
	f=1/y/y
 end function f
!**********************************************************************
 subroutine RK4(x0,xmax,y0,N,xval,yval)
!**********************************************************************
! Negyed rendű Runge-Kutta
! Megoldja a dy/dx=f(x,y) differenciál egyenletet
! az [x0,xmax] tartomanyon, N részre osztva azt.
! A lépésköz h:=(xmax-x0)/(N-1), az eredmény az xval es yval tömb lesz, 
! ami az osztáspontokban számolt x,y értékekket tartalmazza  
!----------------------------------------------------------------------
	implicit none
!  i/o változók
	integer,intent(in):: N
	real(8),intent(in)::x0,xmax,y0
    real(8),intent(out),dimension(:)::xval(N),yval(N)
	real(8)::f
!  Lokális változók
	integer:: i
	real(8)::h,xx,yy,k1,k2,k3,k4
!----------------------------------------------------------------------
   h=(xmax-x0)/(N-1)
   xval(1)=x0;yval(1)=y0
    do i=1,N-1 
     xx = xval(i);yy = yval(i)
     k1 = f(xx    ,yy)                  
     k2 = f(xx+h/2,yy+(h/2)*k1)
     k3 = f(xx+h/2,yy+(h/2)*k2) 
     k4 = f(xx+h  ,yy+h*k3)    
     yval(i+1) = yy + h/6.0d0*(k1+2*k2+2*k3+k4)
	 xval(i+1) = xx+h 
    enddo
   end subroutine RK4