I see a few potential gaps in your methods. One might be that you are using (presumably) constant axes for your aAxis calculation, xAxis and yAxis. I would guess those are the absolute world axes X and Y. That might be ok, if your specific scenario only uses planes oriented to X or Y. That's where you go awry.
To rotate from one direction to another, cross those angles
You cross one plane's normal and an absolute axis. That's not right.
vector3 axis = planeA.normal cross planeB.normal;
However, this axis is only a vector with the same direction as the axis of rotation. A rotation axis needs to be a unit vector.
Cross products need to be normalized
The result of any cross-product, bAxis or other, will only be a unit vector if both crossed vectors are unit vectors AND they are perpindicular to each other. This is a result of this property of cross products (and the fact that sin(90 degrees) = 1):
length(A x B) = sin(angle_between(A, B)) * length(A) * length(B) // so if A and B are unit vectors (have length = 1) length(A x B) = sin(angle_between(A, B)) // and if that length is 1... 1 = sin(angle_between(A, B)) = sin(90 degrees) // A and B were perpindicular
(Side note, you might be able to get a unit vector by crossing a weird combination of non-unit vectors and an angle that coincidentally still adds up to magnitude 1. But that's not the point, and both your planes' normal vectors ought to be unit vectors anyway.)
Rotation matrices are hard to build by hand
If you had normalized your axis, your matrix would still be very wrong. They are not constructed directly from the components of the participating vectors. This is a rotation matrix (again from wikipedia):


Moral of the story: use a library to construct your rotation matrix. You'll thank yourself later.
The correct axis and angle rotation matrix
Good news! You already have it after generating your cross product!
// length(A x B) = sin(angle_between(A, B) // which leaves us with a handy: vector3 a_cross_b = A x B; float angle = arcsign(length(a_cross_b)); vector3 axis = normalize(a_cross_b); matrix correct_rotation_matrix = create_rotation_from_axis_angle(axis, angle); // a quaternion is also acceptable here.
As a bonus, here is some VBA,VB6 code to play with:
'common used data type Type POINT3D X As Double Y As Double Z As Double 'W As Double End Type Function Arcsine(ByVal X As Double) As Double 'inverse sine If X = 1 Then Arcsine = 0.5 * Pi Exit Function End If If X = -1 Then Arcsine = -0.5 * Pi Exit Function End If Arcsine = Sgn(X) * Atn(Sqr(X * X / (1 - X * X))) End Function Function POINT_TRANSFORM(Point As POINT3D, mat() As Double) As POINT3D 'apply transformation matrix to a 3d point Dim P(3) As Double Dim RES As POINT3D P(0) = Point.X P(1) = Point.Y P(2) = Point.Z P(3) = 1 RES.X = P(0) * mat(0, 0) + P(1) * mat(0, 1) + P(2) * mat(0, 2) + P(3) * mat(0, 3) RES.Y = P(0) * mat(1, 0) + P(1) * mat(1, 1) + P(2) * mat(1, 2) + P(3) * mat(1, 3) RES.Z = P(0) * mat(2, 0) + P(1) * mat(2, 1) + P(2) * mat(2, 2) + P(3) * mat(2, 3) 'RES.W = P(0) * mat(3, 0) + P(1) * mat(3, 1) + P(2) * mat(3, 2) + P(3) * mat(3, 3) POINT_TRANSFORM = RES End Function Function VECTOR_CROSS(U As POINT3D, V As POINT3D) As POINT3D 'vector cross product Dim P As POINT3D P.X = ((U.Y * V.Z) - (V.Y * U.Z)) P.Y = ((U.Z * V.X) - (V.Z * U.X)) P.Z = ((U.X * V.Y) - (V.X * U.Y)) VECTOR_CROSS = P End Function Sub CREATE_TM_FROM_AXIS_ANGLE(AXIS As POINT3D, Angle As Double, MATR() As Double) 'axis is given as normalized vector 'angle in RAD 'https://gamedev.stackexchange.com/questions/50880/rotating-plane-to-be-parallel-to-given-normal-via-change-of-basis 'https://en.wikipedia.org/wiki/Rotation_matrix#Axis_and_angle ReDim MATR(3, 3) Dim c As Double Dim S As Double Dim CM As Double X = AXIS.X Y = AXIS.Y Z = AXIS.Z c = Cos(Angle) S = sIn(Angle) CM = 1 - c MATR(0, 0) = X * X * CM + c: MATR(0, 1) = X * Y * CM - Z * S: MATR(0, 2) = X * Z * CM + Y * S MATR(1, 0) = Y * X * CM + Z * S: MATR(1, 1) = Y * Y * CM + c: MATR(1, 2) = Y * Z * CM - X * S MATR(2, 0) = Z * X * CM - Y * S: MATR(2, 1) = Z * Y * CM + X * S: MATR(2, 2) = Z * Z * CM + c MATR(3, 3) = 1 End Sub Function PLANAR_AREAS(A() As Double, B() As Double, MATR() As Double) As Boolean 'given normal vectors aof two planes 'searched: transform matrix to make them coplanar (align them) PLANAR_AREAS = False Dim U As POINT3D Dim V As POINT3D Dim R As POINT3D Dim AXIS As POINT3D Dim L As Double Dim Angle As Double U.X = A(0) U.Y = A(1) U.Z = A(2) V.X = B(0) V.Y = B(1) V.Z = B(2) R = VECTOR_CROSS(U, V) L = Abs(Sqr(R.X * R.X + R.Y * R.Y + R.Z * R.Z)) If L = 0 Then GoTo ulos Angle = 0# * Pi - Arcsin(L) AXIS.X = R.X / L AXIS.Y = R.Y / L AXIS.Z = R.Z / L Call CREATE_TM_FROM_AXIS_ANGLE(AXIS, Angle, MATR()) PLANAR_AREAS = True ulos: End Function Sub TEST_PLANAR_AREAS() 'This routine will make all entitys coplanar to the first ones Dim entity As AcadEntity Dim A() As Double Dim B() As Double Dim MATR() As Double ReDim A(2) ReDim B(2) Dim c As Long Dim V As Variant For Each entity In ThisDrawing.PickfirstSelectionSet 'select a bunch of acad entitys V = entity.Normal 'get their normal vectors If c = 0 Then 'use first ones as reference A(0) = V(0) A(1) = V(1) A(2) = V(2) c = 1 Else 'use rest as "target" B(0) = V(0) B(1) = V(1) B(2) = V(2) If PLANAR_AREAS(A(), B(), MATR()) Then 'calculate transformation matrix entity.TransformBy MATR 'apply Transform Matrix to elements End If End If Next End Sub Sub SHOWNORMALS() 'print the nomals of the autocad entitys (intermediate screen) Dim entity As AcadEntity Dim V Dim I As Long For Each entity In ThisDrawing.PickfirstSelectionSet V = entity.Normal Debug.Print "------" For I = 0 To UBound(V) Debug.Print V(I), Arcsin(CDbl(V(I))) / Pi * 180, Arccos(CDbl(V(I))) / Pi * 180 Next Next End Sub