Ideia:

[m]
A=matrix{4}{3}{
x^2*(1-c)+c x*y*(1-c)-z*sx*z*(1-c)+y*s 0 y*x*(1-c)+z*s y^2*(1-c)+c y*z*(1-c)-x*s 0 x*z*(1-c)-y*s y*z*(1-c)+x*s z^2*(1-c)+c 0 0 0 0 1}
[/m]

[x, y, z] - a vector about which the rotation should be done
c - cos(alpha)

s - sin(alpha)

alpha - rotation angle

Solução:

function CreateGlRotateMatrix(angle, x, y, z: single) : TMatrix;
var
axis: TVector3f;
b, c, ac, s: single;
invLen : Single;
begin

angle:= vectorgeometry.degtorad(angle);

invLen:= RSqrt(x * x + y * y + z * z);
x:= x * invLen;
y:= y * invLen;
z:= z * invLen;

result:= IdentityHmgMatrix;

c:= cos(angle);
s:= sin(angle);

result[0,0] := (x*x) * (1-c)+c;
result[1,0] := x*y * (1-c)-z*s;
result[2,0] := x*z * (1-c)+y*s;

result[0,1] := y*x * (1-c)+z*s;
result[1,1] := (y*y) * (1-c)+c;
result[2,1] := y*z * (1-c)-x*s;

result[0,2] := x*z * (1-c)-y*s;
result[1,2] := y*z * (1-c)+x*s;
result[2,2] := (z*z) * (1-c)+c;
end;

Retirado de: link