c fall with friction eq of motion mdv/dt= -mg +kv**2
      real k,m
      k=8.e-4
      m=0.5
      g=9.8
      vscale=sqrt((m*g)/k)
      tscale=vscale/g
      nstep=2000
      dt=2.*tscale/float(nstep)
c     selects 40 values to be printed
      kprint=int(float(nstep)/40.)
      kount=kprint
c     initial conditions
      x0=0.
      v0=0.
      x1=v0*dt+x0
      v1=dt*(-g+(k/m)*v0**2)+v0
      t=0.
      write(6,100) t, x0,v0
      v0=v1
      do 10 i=2,nstep
      t=dt*float(i)
      v1=dt*(-g+(k/m)*v0**2)+v0
      x2=2.*x1-x0+dt**2*(-g)+(k/m)*(x1-x0)**2
      if(i.eq.kount)then
      kount=kount+kprint
      write(6,100)t,x2,v1
      endif
      x0=x1
      x1=x2
      v0=v1
10    continue
100   format(2x,e10.3,4x,e10.3,4x,e10.3)
      stop
      end
      
      RUN


   0.000E+00     0.000E+00     0.000E+00
   0.399E+00    -0.765E+00    -0.391E+01
   0.799E+00    -0.309E+01    -0.780E+01
   0.120E+01    -0.696E+01    -0.117E+02
   0.160E+01    -0.124E+02    -0.154E+02
   0.200E+01    -0.193E+02    -0.192E+02
   0.240E+01    -0.276E+02    -0.228E+02
   0.280E+01    -0.374E+02    -0.263E+02
   0.319E+01    -0.486E+02    -0.297E+02
   0.359E+01    -0.611E+02    -0.330E+02
   0.399E+01    -0.749E+02    -0.362E+02
   0.439E+01    -0.900E+02    -0.392E+02
   0.479E+01    -0.106E+03    -0.420E+02
   0.519E+01    -0.124E+03    -0.448E+02
   0.559E+01    -0.142E+03    -0.473E+02
   0.599E+01    -0.161E+03    -0.497E+02
   0.639E+01    -0.182E+03    -0.520E+02
   0.679E+01    -0.203E+03    -0.541E+02
   0.719E+01    -0.225E+03    -0.561E+02
   0.759E+01    -0.247E+03    -0.579E+02
   0.799E+01    -0.271E+03    -0.596E+02
   0.839E+01    -0.295E+03    -0.612E+02
   0.878E+01    -0.320E+03    -0.627E+02
   0.918E+01    -0.345E+03    -0.640E+02
   0.958E+01    -0.371E+03    -0.653E+02
   0.998E+01    -0.397E+03    -0.664E+02
   0.104E+02    -0.424E+03    -0.675E+02
   0.108E+02    -0.451E+03    -0.684E+02
   0.112E+02    -0.479E+03    -0.693E+02
   0.116E+02    -0.506E+03    -0.701E+02
   0.120E+02    -0.534E+03    -0.709E+02
   0.124E+02    -0.563E+03    -0.715E+02
   0.128E+02    -0.592E+03    -0.721E+02
   0.132E+02    -0.621E+03    -0.727E+02
   0.136E+02    -0.650E+03    -0.732E+02
   0.140E+02    -0.679E+03    -0.737E+02
   0.144E+02    -0.708E+03    -0.741E+02
   0.148E+02    -0.738E+03    -0.745E+02
   0.152E+02    -0.768E+03    -0.748E+02
   0.156E+02    -0.798E+03    -0.752E+02
   0.160E+02    -0.828E+03    -0.755E+02