c  test of runge-kutta  method
      real t,x(4),dt,t0,x0(4)
      real xreal,error
      real s,s2
      real tlim
      real pi
      integer n,nwrite, i
      common /para/alpha
c
      open(unit=2,file='xy.dat')
      pi=3.1415926535897932385
      write(6,*) ' v0, theta ?'
      read(5,*) v0,theta
      write(6,*) ' alpha ?'
      read(5,*) alpha
      t0=0.0                                                            
      x0(1)=0.0
      x0(2)=0.0
      x0(3)=v0*cos(pi*theta/180)
      x0(4)=v0*sin(pi*theta/180)
      write(6,*) ' dt ?'
      read(5,*) dt
      i=0
      t=t0                                                              
      x(1)=x0(1)                        
      x(2)=x0(2)                        
      x(3)=x0(3)                        
      x(4)=x0(4)                        
   10 continue
      i=i+1
      call runge2(4,t,x,dt)                                               
      write(2,601) t,x(1),x(2),x(3),x(4)                                     
      write(6,601) t,x(1),x(2),x(3),x(4)                                     
      if(x(2).lt.0.0) then
        stop
      endif
      go to 10                                                          
  601 format(1h ,5(1pe12.5,2x))                                         
      end                                                               
*********************************************************************** 
      real function f(i,t,x)                                                   
      integer i
      real t
      real x(4)
      real pi
      common /para/alpha
      pi=3.1415926535897932385
      if(i.eq.1) then	
        f=x(3)                                          
      elseif(i.eq.2) then
        f=x(4)
      elseif(i.eq.3) then
        f=0.0 -alpha*x(3)
      elseif(i.eq.4) then
        f=-9.8 -alpha*x(4)
      endif
      return                                                            
      end                                                               
*********************************************************************** 
c  subroutine of runge-kutta method of many valuables
      subroutine runge2(n,x,y,dx)                                       
      real h1,k(10,4),k1(10),y1(10),y2(10),y3(10),y(10)
c
      h1=dx/2.0                                                         
      do 10 i1=1,n
        k(i1,1)=f(i1,x,y)*dx                                              
        y1(i1)=y(i1)+k(i1,1)/2.0
  10  continue                                                          
      do 20 i2=1,n
        k(i2,2)=f(i2,x+h1,y1)*dx                                          
        y2(i2)=y(i2)+k(i2,2)/2.0
  20  continue                                                          
      do 30 i3=1,n
        k(i3,3)=f(i3,x+h1,y2)*dx                                          
        y3(i3)=y(i3)+k(i3,3)
  30  continue                                                          
      do 40 i4=1,n
        k(i4,4)=f(i4,x+dx,y3)*dx                                          
   40 continue                                      
c
      x=x+dx                                                            
      do 50 i5=1,n
        k1(i5)=(k(i5,1)+2.0*k(i5,2)+2.0*k(i5,3)+k(i5,4))/6.0              
        y(i5)=y(i5)+k1(i5)
  50  continue                                                          
      return
      end