A = [-1 0 8; 0 -2 0; 0 0 -3]
B = [0 2; 3 0; 0 0]
C = [1 0 -1]
D = [0 0]

%

[M, L] = eig(A)

%

syms t
Phi_p = @(t) diag(exp(diag(L) * t));
Phi_p(t)

%

M_inv = inv(M);
Phi = @(t) M * Phi_p(t) * M_inv;
Phi(t)

%

u = @(t) [4; sin(2 * pi * t)];
u(t)

%

syms T
x_fo = @(t) int(Phi(t - T) * B * u(T), T, 0, t);
x_fo(t)

%

y_fo = @(t) C * x_fo(t) + D * u(t);
y_fo(t)

%

t_vec = 0:0.02:7;
y_fun = matlabFunction(y_fo(t));
plot(t_vec, y_fun(t_vec))