So, I may have recreated the wheel, but I couldn’t find the code for it online. 
This is a function which uses two objects as parameters, finds their locations and then rotates
the Z-axis, then Y-axis of the rot_obj parameter accordingly so that the Z-axis of rot_obj is pointing straight at the selected object. THIS IS ALMOST EXACTLY LIKE THE TRACK-TO CONSTRAINT IN BGE.
IS THERE A BETTER WAY? I felt like pulling my hair out while creating this. lol
def faceTowards(selected,rot_obj): #Each location is a vector.
vec1 = mathutils.Vector((selected.location))
vec2 = mathutils.Vector((rot_obj.location))
#print("Vec1",vec1,"Vec2", vec2)
#print(findDistance(vec2,vec1))
matrix = rot_obj.matrix_world #RESETS SCALE AND ROTATION
matrix[0]=(1,0,0,0)
matrix[1]=(0,1,0,0)
matrix[2]=(0,0,1,0)
matrix[3]=(0,0,0,1)
newvector = vec1-vec2
#Delta -> aka: DIFFERENCE or 'Distance' in this case
deltax = newvector[0] #Distance between two x locations... etc
deltay = newvector[1]
deltaz = newvector[2]
#ATAN or ArcTan = adjacent/opposite
#The following 'angle' variables are the angles between associated Distances.
#Think of the distances(delta(x,y,z)) as now being lines with a length of each distance.
angle1 = math.atan2(deltay,deltaz)#The resulting angle variable is in radians
angle2 = math.atan2(deltax,deltaz)
angle3 = math.atan2(deltay,deltax)
new_rotation = [] #The array we turn into a vector with the new rotations.
#The following works but was found via trial and error results, so I can't explain very well.
new_rotation.append(((angle1)+math.radians(180))) #IGNORE THIS X AXIS AS WE ARE USING Y AXIS ROTATION.
measure = ((angle2)+math.radians(180))
if measure>=math.radians(180):
measure-=math.radians(360)
measure = -measure
new_rotation.append(measure+math.radians(180)) ###THIS IS FOR Y AXIS ROTATION
new_rotation.append(((angle3)+math.radians(180))) ###THIS IS FOR Z AXIS ROTATION
###LITTLE PROGRAM FOR DISPLAYING THE ANGLES IN THIS ORDER: FIRST, RADIANS, THEN, DEGREES
#count = 0
#for x in new_rotation:
#print("Axis: ",count)
#count+=1
# print(x,math.degrees(x))
########################################################################
ROTATION = matrix.to_euler('XYZ')
ROTATION.rotate_axis('Z',new_rotation[2]) ## Z-AXIS FACES POINT... NOTE: Z-AXIS IS 'NORMAL' HERE
#test.rotate_axis('X',faceTowards(vec2,vec1)[0])
ROTATION.rotate_axis('Y',new_rotation[1]) ###LOCAL Y-AXIS
ROTATION =ROTATION.to_matrix().to_4x4()
rot_obj.matrix_world=ROTATION
return print("Rotation: Done")
testroom.blend (909 KB)