c
c 3414 Summer 2002
c HW4, Problem 2(c)
c

c
c taken from page 315-316: numerical mathematics and computing, 
c  cheney/kincaid, 1985
c
c file: rk4.for
c
c runge-kutta method of order 4 for solving an initial value problem (rk4,f)
c
      external f
c     data t/0.0/, x/16.0/, h/0.1/, nstep/20/
      data t/0.0/, x/16.0/
      read *, nstep
      h = 2.0/nstep
      call rk4(f,t,x,h,nstep) 
      stop
      end 
  
      function f(t,v) 
      real   m, g, k

      m = 0.11
      g = -9.8
      k = 0.002

      f = g - k*v*abs(v)/m
      return
      end
  
      subroutine rk4(f,t,x,h,nstep)   
      print 3,t,x 
      h2 = 0.5*h
      start = t   
      do 2 k = 1,nstep      
        f1 = h*f(t,x)       
        f2 = h*f(t + h2,x + 0.5*f1)   
        f3 = h*f(t + h2,x + 0.5*f2)   
        f4 = h*f(t + h,x + f3)
        x = x + (f1 + f2 + f2  + f3 + f3 + f4)/6.0
        t = start + h*real(k) 
        print 3,t,x 
   2  continue
   3  format(f5.2,f8.4)
      return
      end