import pylab

delta_t = 1e-3 # 1 ms

t = 0.0
b = 6.0
M = 25.0
f = 30 # newton

v = 0.0
p = 0.0

vettore_pos = [ ]
vettore_vel = [ ]
vettore_tempi = [ ]

while t < 20:
    p = p + delta_t * v
    v = (1 - b*delta_t/M) * v + delta_t / M * f
    t = t + delta_t
    vettore_vel.append(v)
    vettore_pos.append(p)
    vettore_tempi.append(t)


pylab.figure(1)
pylab.plot(vettore_tempi, vettore_vel, 'r-+', label='vel, v(t)')
pylab.xlabel('time')
pylab.legend()

pylab.figure(2)
pylab.plot(vettore_tempi, vettore_pos, 'b-+', label='position, p(t)')
pylab.xlabel('time')
pylab.legend()

pylab.show()