c  periodic force
       realk,m
       equivalence (dt,h)
       ampl=2.
       period=10.
       k=1.
       m=1.
       b=0.2
       omega=sqrt(k/m)
       nstep=2000
       t1=1./b
       t2=1./omega
       t3=period
       tscale=min(t1,t2,t3)
       dt=20.*tscale/float(nstep)
      x0=1.
      v0=0.
      x1=v0*dt+x0
      kprint=int(float(nstep)/100.)
      kount=kprint
      t=0.
      write(6,100)t,x0,fuerza(0.,period,ampl)
      do 10 i=2,nstep
      t=float(i)*h
      x2=2.*x1-x0-dt*(b/m)*(x1-x0)+dt**2*(-(k/m)*x1+
     $fuerza(t-dt,period,ampl)/m)
      v=(x2-x1)/dt
      if(i.eq.kount)then
      write(6,100)t,x2,fuerza(t,period,ampl)
      kount=kount+kprint
      endif
      x0=x1
      x1=x2
10    continue
100   format(3x,'t,x,force=',3x,e10.3,3x,e10.3,3x,e10.3)
      stop
      end
      
      function fuerza(t,period,ampl)
      half=period/2.
      ifact=int(t/period)
      tprime=t-float(ifact)*period
      if(tprime.le.half)fuerza=+ampl
      if(tprime.gt.half)fuerza= -ampl
      return
      end