J = 1;
K = 9;
b = 0.6;

%

A = [-b / J, 1 / J; -K, 0]
B = [0; K]

%

[M, D] = eig(A)

%

dt = 0.2;
Phi_p = diag(exp(diag(D * dt)))
Phi = real(M * Phi_p * inv(M))

%

real(M * diag(exp(diag(D * 2 * dt))) * inv(M))

%

Phi^2

%

syms t T
u = 50;
M_inv = inv(M);
Phi = @(t) M * diag(exp(diag(D * t))) * M_inv;
x_fo = @(t) int(Phi(t - T) * B * u, T, 0, t);
x_fo(t)

%

t_vec = 0:0.05:15;
x_fun = matlabFunction(x_fo(t));
plot(t_vec, x_fun(t_vec))

%

syms tu
u = @(t) subs(piecewise(tu<10, 5 * tu, tu>=10, 50), tu, t);
x_fo = @(t) int(Phi(t - T) * B * u(T), T, 0, t);
x_fo(t)

%

x_fun = @(t) subs(x_fo(tu), tu, t)

%

plot(t_vec, x_fun(t_vec))
hold on
plot(t_vec, u(t_vec))