Ok i solveI have solved the problem myself. In summary, one has to evaluate the complex matrix for each argument of the Fourier/Laplace transform $u_1, \dots, u_{1000}$. This cannot be done with MatLab's element-by-element product. So if anyoneI used a for loop that would multiply for each argument the matrix by the scalar $u_1, \dots, u_{1000}$.
Here is interested, i leave the correctrevised code:
function resu = chf_wmsv(u,t,r,q,M,R,Q,Sigma0,S0,beta)
i = complex(0,1);
resu = zeros(1,length(u));
for idx = 1:length(u)
% Exponential of the matrix (14) in Da Fonseca et al. [2008]
w = u(1,idx);
MATEXP = expm(t*[M, -2.0*(Q'*Q); ...
0.5*i*w*(i*w - 1)*eye(2), -1.0*(M' + 2.0*i*w*R*Q)]);
% A11 = [MATEXP(1:2,1:2)];
% A12 = [MATEXP(1:2,3:4)];
A21 = [MATEXP(3:4,1:2)];
A22 = [MATEXP(3:4,3:4)];
% Computation of matrix function A(tau)
A = A22\A21;
% Computation of scalar function C(tau)
C = -0.5*beta * trace(log(A22) + t*(M' + 2*R*Q)) + i*w*(r-q);
% Characteristic Function for WMSV model
resu(1,idx) = exp(trace(A*Sigma0) + C + i*w*log(S0));
end
end
