Post your Mult function.

EDIT:

My guess is you've made a mistake in Mult, also possibly not normalizing in the right location.

------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------

The code below is the Quaternion implementation from Octave 3.2.4 (basically open source MatLab). The syntax is simple so hopefully the language difference won't be a problem as this is a math question. I removed the error checking portions from the code but attached the full source in a zip file for you.

Code:

function v = qtrans (v, q)
v = qmult (q, qmult (v, qinv (q)));
endfunction

Code:

function retval = qmult (a, b)
[a1, b1, c1, d1] = quaternion (a);
[a2, b2, c2, d2] = quaternion (b);
ri = b1*c2 - c1*b2 + d1*a2 + a1*d2;
rj = c1*a2 - a1*c2 + d1*b2 + b1*d2;
rk = a1*b2 - b1*a2 + d1*c2 + c1*d2;
rr = -(a1*a2 + b1*b2 + c1*c2) + d1*d2;
retval = quaternion (ri, rj, rk, rr);
endfunction

Code:

function retval = qinv (q)
if (norm (q) != 0)
retval = qconj (q) / sum (q .* q);
else
error ("qinv: zero quaternion passed!");
endif
endfunction

Code:

function retval = qconj (q)
[a, b, c, d] = quaternion (q);
retval = quaternion (-a, -b, -c, d);
endfunction