Tough one, eh?  :-)  I've googled all over the place, and I've found some code, but it does not work. that is, the resulting quaternion, when used as obj.rotate_quaternion, does not result in properly rotated cylinder with ep1 = (0,0,0) and ep1 = (0,100,0)

def FindRotation(ep1, ep2):

First, create quaternion per

http://stackoverflow.com/questions/1171849/finding-quaternion-representing-the-rotation-from-one-vector-to-another

vec = numpy.array(ep2) - numpy.array(ep1)

vec = vec/numpy.linalg.norm(vec)

if not vec[0] and not vec[1]:
return (0,0,0)
print ("vec: %s"%str(vec))
Z = numpy.array((0,0,1))
kcostheta = numpy.dot(Z, vec)
k = Length(vec) # sqrt(Length(vec)^2 * Length(Z)^2) = Length(vec)
Q = Quaternion(numpy.cross(vec,Z), kcostheta+k)
Q.normalize()
print ("kcostheta: %f, k: %f, Q: %s"%(kcostheta, k, str(Q)))
return Q