How can I create a General force integration

I am trying to create a program for a multi degree freedom code. I am using the modal matrix approach and for that I must integrate the general force for a vibration on a system. I have wrote it down below.

The full equation which I am trying to find is:

The problem is I am not sure how to incorporate Q(tao) in my integral. For this particular problem it is a constant. I think matlab thinks of tao as a row instead of symbol. Also tao is p in the code. Thanks

Here is my code:

%Code for multi-degree of freedom
clear all

%input mass matrix
m = 1; %we can change this accordingly 

M = [1,0,0;0,1,0;0,0,1] * m;

M_inv= inv(M);

%input stiffness matrix
k= 1; %we can also change this

K = [2,-1,0;-1,2,-1;0,-1,1] *k;

%Now for eigenvalue problem
MK = M\K;

[V,D] = eig(MK);

[eig_MK,ind] = sort(diag(D));
% Ds = D(ind,ind);
% Vs = V(:,ind);

w1 = sqrt(eig_MK(1));
w2 = sqrt(eig_MK(2));
w3 = sqrt(eig_MK(3));

t=0;
F1i=0;
F2i=0;
F3i=0;

wf1=0;
wf2=0;
wf3=0;

%Initial conditions
t=0;
x0= [1;0;0];
v0 = [0;0;0];

qx0= V'* M * x0;
qv0 = V'* M * v0;


F1 = F1i*cos(wf1*t);
F2 = F2i*cos(wf2*t);
F3 = F3i*cos(wf3*t);

Ft= [F1;F2;F3];

Qt= V'*Ft;
i=1;
syms p
for t=1:100
% Have to mae Qf2 and Qf3 and qt2 and qt3
Qf1(i) = Qt(1);
I= (1/w1)*int(Qf1(p)*sin(w1*(t-p)), p, 0, t);
qt1(i) = qx0(1)*cos(w1*t) + (qv0(1)/w1)*sin(w1*t)+I;
i=i+1;

end