% compute the tf, the char. poly. and poles of a LTI system,
% in particular the mass-spring-damper system, 
clear all
clc
close all

% define the symbolic variables
syms s k b m

% define A, B, C, D in parametric form
A=[0 1; -k/m -b/m ];     
B=[0;1/m]; 
C=[ 1 0]; 
D=0; 

I=eye(2);

% compute the transition matrix 
Phi_s=inv(s*I-A);

% compute the tf
G_s=C*Phi_s*B+D;

% compute char poly.
D_s=charpoly(A);
% determine the coeff. of the char. equation s^2+a1*s+a0=0
D_s_n=D_s/D_s(1);
a0=D_s_n(end);
a1=D_s_n(end-1);

% determine zita and omega_n
omega_n=sqrt(a0);
zita=a1/2/sqrt(a0);

% determine b to get zita =1
zita_eq_one=zita-1;
b_zita_eq_one=solve(zita_eq_one,b);

% value for k, k=4
b1_val_m=subs(b_zita_eq_one,k,4);
% value for  m, m=1
b1_val=subs(b1_val_m,m,1);

% or
D_s_1=det(s*I-A);

% compute poles
poles_G=roots(D_s);

% or define char. eq. and compute the roots
% define char. eq.
eq_c=D_s_1==0;
% then compute the roots, i.e. poles
poles_eqc=solve(eq_c,'ReturnConditions',true);
poles_sol=poles_eqc.s;


% compute Delta of char. eq
delta_eq=D_s(2)^2-4*D_s(1)*D_s(3);

% find parameters (in particolar b) to get Delta=0;
delta_zero=solve(delta_eq,b);

% value for k, k=4
b2_val_m=subs(delta_zero,k,4);

% value for  m, m=1
b2_val=subs(b2_val_m,m,1);