c damped oscillations -see Jerry Marion-Classical Dynamics- p.121
real k,m
k=1.
m=1.
A=1.
beta=0.2
omega=sqrt(k/m)
nstep=2000
t1=1./beta
t2=1./omega
tscale=min(t1,t2)
dt=13.*tscale/float(nstep)
c initial conditions
x0=0.
v0=a*sqrt(omega**2-beta**2)
x1=v0*dt+x0
kprint=int(float(nstep)/50.)
kount=kprint
t=0.
write(6,100)t,x0,v0
do 10 i=2,nstep
t=dt*float(i)
vel=(x1-x0)/dt
x2=2.*x1-x0+dt**2*(-beta*vel-omega**2*x1)
if(i.eq.kount)then
v=(x2-x1)/dt
write(6,100)t,x2,v
kount=kount+kprint
endif
x0=x1
x1=x2
10 continue
100 format(2x,'t,x,v=',3x,e10.3,3x,e10.3,3x,e10.3)
stop
end
RUN
t,x,v= 0.000E+00 0.000E+00 0.980E+00
t,x,v= 0.260E+00 0.246E+00 0.900E+00
t,x,v= 0.520E+00 0.463E+00 0.765E+00
t,x,v= 0.780E+00 0.638E+00 0.586E+00
t,x,v= 0.104E+01 0.763E+00 0.377E+00
t,x,v= 0.130E+01 0.832E+00 0.155E+00
t,x,v= 0.156E+01 0.843E+00 -0.664E-01
t,x,v= 0.182E+01 0.797E+00 -0.272E+00
t,x,v= 0.208E+01 0.702E+00 -0.450E+00
t,x,v= 0.234E+01 0.566E+00 -0.589E+00
t,x,v= 0.260E+01 0.399E+00 -0.682E+00
t,x,v= 0.286E+01 0.215E+00 -0.726E+00
t,x,v= 0.312E+01 0.261E-01 -0.720E+00
t,x,v= 0.338E+01 -0.155E+00 -0.667E+00
t,x,v= 0.364E+01 -0.317E+00 -0.573E+00
t,x,v= 0.390E+01 -0.450E+00 -0.447E+00
t,x,v= 0.416E+01 -0.546E+00 -0.297E+00
t,x,v= 0.442E+01 -0.602E+00 -0.136E+00
t,x,v= 0.468E+01 -0.616E+00 0.263E-01
t,x,v= 0.494E+01 -0.588E+00 0.178E+00
t,x,v= 0.520E+01 -0.524E+00 0.311E+00
t,x,v= 0.546E+01 -0.428E+00 0.417E+00
t,x,v= 0.572E+01 -0.309E+00 0.490E+00
t,x,v= 0.598E+01 -0.176E+00 0.527E+00
t,x,v= 0.624E+01 -0.382E-01 0.528E+00
t,x,v= 0.650E+01 0.952E-01 0.494E+00
t,x,v= 0.676E+01 0.216E+00 0.429E+00
t,x,v= 0.702E+01 0.316E+00 0.340E+00
t,x,v= 0.728E+01 0.390E+00 0.233E+00
t,x,v= 0.754E+01 0.435E+00 0.116E+00
t,x,v= 0.780E+01 0.449E+00 -0.287E-02
t,x,v= 0.806E+01 0.434E+00 -0.115E+00
t,x,v= 0.832E+01 0.390E+00 -0.214E+00
t,x,v= 0.858E+01 0.323E+00 -0.295E+00
t,x,v= 0.884E+01 0.238E+00 -0.351E+00
t,x,v= 0.910E+01 0.142E+00 -0.382E+00
t,x,v= 0.936E+01 0.419E-01 -0.386E+00
t,x,v= 0.962E+01 -0.562E-01 -0.365E+00
t,x,v= 0.988E+01 -0.146E+00 -0.321E+00
t,x,v= 0.101E+02 -0.221E+00 -0.258E+00
t,x,v= 0.104E+02 -0.278E+00 -0.181E+00
t,x,v= 0.107E+02 -0.314E+00 -0.966E-01
t,x,v= 0.109E+02 -0.328E+00 -0.987E-02
t,x,v= 0.112E+02 -0.319E+00 0.731E-01
t,x,v= 0.114E+02 -0.290E+00 0.147E+00
t,x,v= 0.117E+02 -0.243E+00 0.208E+00
t,x,v= 0.120E+02 -0.183E+00 0.251E+00
t,x,v= 0.122E+02 -0.114E+00 0.277E+00
t,x,v= 0.125E+02 -0.408E-01 0.282E+00
t,x,v= 0.127E+02 0.312E-01 0.269E+00
t,x,v= 0.130E+02 0.976E-01 0.239E+00