Aligning plane normal vector side to face a point : Python
The function precisely for that is to_track_quat 
Works much like track to constraints. Gives more than rotation difference.

to_track_quat(track, up)


This uses the Vector to point with the specified track axis, (Z for normal) and a secondary "up" axis (x or y, never same ...).

Code would be:

def align_plane_to_point(obj, point):
    dir = point - obj.location
    return dir.to_track_quat("Z", "X").to_euler()


if obj:
    obj.rotation_euler = align_plane_to_point(obj, point)


Of course you can extra check for obj rotation type, the function giving a quaternion, I convert it here to euler assuming the usual kind. Or may use it for rotation matrix, like:

rotMatrix = dir.to_track_quat("Z", "X").to_matrix().to_4x4()


(as quat to matrix produces a 3x3 matrix)


If u need more control, but more complicated, you can have the second dir (up be controlled by another vector) you may look into this, if it helps.


I build there the rotation matrix directly out of 1 vector and a second that is not necessarily perpendicular. 

First I make the second/third vector perpendicular by the cross functions, with some exceptions. I'll not copy that here ... Then can build a 3x3 matrix (only rotation) by putting the 3 vecs into columns

mat3x3 = mathutils.Matrix().to_3x3()"
mat3x3.col[0], mat3x3.col[1], mat3x3.col[2] = mx, my, mz"
eulerRotation = mat3x3.to_euler()


where 3 vecs,  mx, my, mz,  mz would be the point - obj.location, the other 2 being determined as in the linked code

if u use directly such matrix (to 4x4 ok) you'll have some difficulty controlling the scale for some exceptional cases (may fall to 0 on some axes), so to euler or to quaternion is the way, even if u further convert back to a full matrix