It seem like the script only apply a transform matrix to object origin, which doesn't consider the local location of the face center.
You will need to calculate the difference vector between that face center and the origin of the same object, and move that object by the vector.
Here is the script I modified (with mess, Blender 2.8 API):
import bpy
from mathutils import Matrix, Vector
import bmesh
bpy.ops.object.mode_set(mode='EDIT')
context = bpy.context
scene = context.scene
A = context.object
A.select_set(state=False)
B = bpy.context.selected_objects[0]
A.select_set(state=True)
#obj = context.edit_object (!) equal to A
src_mw = A.matrix_world.copy()
src_bm = bmesh.from_edit_mesh(A.data)
src_face = src_bm.select_history.active
src_o = src_face.calc_center_median()
src_normal = src_face.normal
src_tan = src_face.calc_tangent_edge()
# This is the target, we change the sign of normal to stick face to face
dst_mw = B.matrix_world.copy()
dst_bm = bmesh.from_edit_mesh(B.data)
dst_face = dst_bm.select_history.active
dst_o = dst_face.calc_center_median()
dst_normal = -(dst_face.normal)
dst_tan = (dst_face.calc_tangent_edge())
vec2 = src_normal @ src_mw.inverted()
matrix_rotate = dst_normal.rotation_difference(vec2).to_matrix().to_4x4()
vec1 = src_tan @ src_mw.inverted()
dst_tan = dst_tan @ matrix_rotate.inverted()
mat_tmp = dst_tan.rotation_difference(vec1).to_matrix().to_4x4()
matrix_rotate = mat_tmp @ matrix_rotate
matrix_translation = Matrix.Translation(src_mw @ src_o)
# This line applied the matrix_translation and matrix_rotate
B.matrix_world = matrix_translation @ matrix_rotate.to_4x4()
# We need to recalculate these value since we change the matrix_world
dst_mw = B.matrix_world.copy()
dst_bm = bmesh.from_edit_mesh(B.data)
dst_face = dst_bm.select_history.active
dst_o = dst_face.calc_center_median()
# The following is telling blender to find a translation from face center to origin,
# And than apply it on world matrix
# Be Careful, the order of the matrix multiplication change the result,
# We always put the transform matrix on "Left Hand Side" to perform the task
dif_mat = Matrix.Translation(B.location - dst_mw @ dst_o)
B.matrix_world = dif_mat @ B.matrix_world
And since the face in quad is actually not flat in most of the case, you may find out the script is not always stick the face beautifully because it is impossible to bend that face with fixed vertices position.
You may want to split the polygon in to triangle to make the face flat.
And for the matrix tutorials, it's just math problem in 3-dimension space. You may want to check out some guide of how these surface model work in most of the 3D application.
Matrix Transformation: 3D Transformation - tutorialspoint