I would use a pesudoinverse:
A = [6 4;6 1; 1 2; 6 4]
B = [8 8 6 7 5; 8 8 1 6 0; 1 4 3 8 7]
C = [18042 21288 10716 22446 12924; 12768 15024 7593 15795 9099;...
5351 6332 3174 6697 3854; 18042 21288 10716 22446 12924]
How about thinking of Moore-Penrose pseudoinverse?
X = pinv(A)*C*pinv(B)
Check if this works well:
A*X*B - C
Floating-point relative accuracy is given by
eps