Atan2 works great for 2D but the question is in 3D.
3D requires up to 3 angles to be computed.
I have handled this a couple different ways depending on my needs at the time.
The 'simplest' is to use Atan2 but you need to calculate two angles (One along the horizontal 'ground' plane, and one along the vertical)
So, assuming you are using xna-space (z-x is the ground plane, y is 'up') you can do something like this
Vector3 A,B,deltaAB;//set A,B how ever you need
deltaAB=B-A;
float horizontalAngle=(float)Math.ATan2(deltaAB.X,deltaAB.Z);
//to calculate the vertical angle use the length of the ground component of deltaAB and deltaAB.Z
Vector3 deltaAB_Ground=deltaAB;
deltaAB_Ground=0;//make flat
float verticalAngle=(float)Math.ATan2(deltaAB.Z,deltaAB_Ground.Length());
one way to use these angles (depending on how you want to define your rotations) is
Matrix horTransform=Matrix.CreateRotationY(horizontalAngle);
Matrix verTransform=Matrix.CreateRotatoinZ(verticalAngle);//this could be CreateRotationX instead depending on which 'vertical' rotatoin you want.
the other method flows simialr to this (and again it depends on your particular needs):
Vector3 A,B,deltaAB,unitAB;
deltaAB=B-A;
unitAB=Vector3.Normalize(deltaAB);
Matrix transform=Matrix.CreateWorld(new Vector3(0,0,0),unitAB,new Vector3(0,1,0));//this will give you a single matrix that applies a rotation on all 3-axis based on the direction of unitAB. Also allows for a translation and definiation of 'up'
code with pride, not ego!