c harmonic oscillator (classical)- Fisi 3011 august 2005
      real k ,mass
      a(x,v)=-(k/mass)*x -(b/mass)*v
      data k,mass,b /1.,1.,0./
      data x0,v/1.,0./
      tscale=sqrt(mass/k)
      nstep=3000
      dt=6.5*tscale/float(nstep)
      x1=v*dt+x0
c kprint and kount defines print interval
      kprint=int(float(nstep)/80.)
      kount=kprint
      t=0.
      print 100,t,x0
      do 10 i=2,nstep
      t=dt*float(i)
      x2=2.*x1-x0+dt**2*a(x1,v)
      v=(x2-x1)/dt
      if(i.eq.kount)then
      print 100 ,t,x2
      kount=kount+kprint
      endif
      x0=x1
      x1=x2
10    continue
100   format(2x,'t, x=',2(4x,e11.4))
      stop
      end