Convolution integral with matrices

Hey everybody! My question is quite simple I guess:

I would like to perform the following convolution integral (written in latex form for the sake of clarity):

, where A is 2x2, B is 2x1 and u is scalar.

I would like to compute this WITHOUT symbolic functions.

So far, this is what I've done:

T = 0:0.01:20 % Assume 20 seconds of simulation.

% Compute the values of expm(A) and B*u for each time step ---------
  for i=1:length(T)
    tt=t(i);
    E(:,:,i) = expm(A*tt);
    U(:,i) = B*u(tt);
  end

 nu = size(A,2); % Assume that A is 2x2.
 result2 = 0;
 for j=1:nu
    AA=squeeze(A(:,j,:));
    result2 = result2 + median(diff(T))*conv2(AA,B(j,:),'full'); 
 end

This obtained result is not correct, and I'm struggling to find why.

Thanks!

[Niraj Choudhary]

@Nicolás Faedo have u got the answer?