Usually the best way to think about this is as rotation about an axis. Consider this example:
origin = [0 0 0];
pt1 = randn(1,3);
pt2 = randn(1,3);
scatter3(origin(1),origin(2),origin(3),'Marker','o')
hold on
line([origin(1) pt1(1)],[origin(2) pt1(2)],[origin(3) pt1(3)],'Color','red')
line([origin(1) pt2(1)],[origin(2) pt2(2)],[origin(3) pt2(3)],'Color','green')
axis vis3d
axis equal
xlim([-2 2])
ylim([-2 2])
zlim([-2 2])
We want an axis which is normal to the plane which passes through our two points and the origin. We can get that with a cross product:
u = pt1-origin;
v = pt2-origin;
u = u/norm(u);
v = v/norm(v);
n = cross(u,v);
n = n/norm(n)
line([origin(1) n(1)],[origin(2) n(2)],[origin(3) n(3)],'Color',[.75 .75 .75])
Now we can generate rotation matrices like so:
makehgtform('axisrotate',n,ang);
So all we need is the correct value for ang. The dot product of the two vectors is the cosine of the angle:
g = hgtransform;
line([origin(1) pt1(1)], ...
[origin(2) pt1(2)], ...
[origin(3) pt1(3)], ...
'Color','blue','LineStyle',':','LineWidth',2,'Parent',g)
ang = acos(dot(u,v));
g.Matrix = makehgtform('axisrotate',n,ang);
The one thing to watch out for he is that you need to do a little bookkeeping with the sign of your dot product. Sometimes you'll find yourself going around the axis the "wrong" way.
Does that make sense?