Instead of using:

trans_mat := unapply(<rot_mat(<x[1], x[2], x[3]>, theta) | (LinearAlgebra[IdentityMatrix](3) - rot_mat(<x[1], x[2], x[3]>, theta)) . <p[1], p[2], p[3]>>, x::Vector, p::Vector, theta)

, use

trans_mat := (x::Vector, p::Vector, theta) -> <rot_mat(<x[1], x[2], x[3]>, theta) | (LinearAlgebra[IdentityMatrix](3) - rot_mat(<x[1], x[2], x[3]>, theta)) . <p[1], p[2], p[3]>>

This creates a proper matrix. On further examination of the output from my original post, it appears thatis just a concatenation of elements, with no regards to making it into a combined matrix of 3x4 size.

This answers my original question, which is

trans_mat_ang := (x::Vector, p::Vector, theta) -> <<rot_mat(<x[1], x[2], x[3]>, theta) | (LinearAlgebra[IdentityMatrix](3) - rot_mat(<x[1], x[2], x[3]>, theta)) . <p[1], p[2], p[3]>>, <0 | 0 | 0 | 1>>;

@Rouben Rostamian  Thank you very much for your reply. Indeed, from a functions' point of view, your suggestion is more elegant. But the reason I have to stick with Homogeneous transforms is because it makes it easier to deal with appyling Screw Theory to Robotics. I learned it from this book, and since it was already making life easier to deal with the dynamics of multi-rigid-body systems, I will have to stick with it till I get more comfortable without it.