c  test of runge-kutta  method
      real x,y,dx,x0,y0
      real yreal,error
      real s,s2
      real xlim
      real pi
      integer n,nwrite, i
c
      pi=3.1415926535897932385
      x0=0.0                                                            
      y0=0.0
      xlim=1.0
      write(6,*) ' dx ?'
      read(5,*) dx
      n=nint(xlim/dx)
      nwrite=n/10
      if(nwrite.eq.0) nwrite=1
      i=0
      s2=0.0
      x=x0                                                              
      y=y0                        
   10 continue
      i=i+1
      call runge1(x,y,dx)                                               
      yreal=sin(2.0*pi*x)
      error=y-yreal
      s2=s2+error**2
      write(2,601) x,y,yreal,error                                             
      if(mod(i,nwrite).eq.0) then
        write(6,601) x,y,yreal,error                                           
      endif
      if(x.gt.xlim) then
        s2=s2/i
        s=sqrt(s2)                                                 
        write(6,*) ' Number of Steps= ',i,'  sigma= ',s
        stop
      endif
      go to 10                                                          
  601 format(1h ,4(1pe12.5,2x))                                         
      end                                                               
*********************************************************************** 
      real function f(x,y)                                                   
      real x,y
      real pi
      pi=3.1415926535897932385
      f=2.0*pi*cos(2.0*3.1415926*x)                                           
      return                                                            
      end                                                               
*********************************************************************** 
c  subroutine of runge -kutta method of one variable
      subroutine runge1(x,y,dx)                                         
      real x,y,dx
      real  h1,k(4),k1
c
      h1=dx/2.0                                                         
      k(1)=f(x,y)*dx
      y1=y+k(1)/2.0                                                     
      k(2)=f(x+h1,y1)*dx
      y2=y+k(2)/2.0                                                     
      k(3)=f(x+h1,y2)*dx
      y3=y+k(3)                                                         
      k(4)=f(x+dx,y3)*dx
c
      x=x+dx                                                            
      k1=(k(1)+2.0*k(2)+2.0*k(3)+k(4))/6.0
      y=y+k1                                                            
      return
      end