Solving this matrix with the correct response using MATLAB and simulink

this is the given problem.

this is the Simulink and here is the call back for the function command:

here is my code

clc

clear

%---------------EQUATIONS_FROM_FBD------------------------------------------------------------------------------------------------------------------------------------------------

%......................................................................SUM_OF_EXTERNAL_FORCES=M*X"................................................................................

% MfXf"=-Ksf(Xf-[Xb-(L1*theta)})-Bsf(Xf'-{Xb'-(L1*theta')})+(Kf*Xf);

% Vf'=(1/Mf)[(-Ksf(Xf-(L1*theta))-(Bsf(Vf-(L1*theta'))+(Kf*Xf)] ;

% MrXr"=-Ksr(Xr+[Xb+(L2*theta)])-Bsr(Xr'+[Xb'+(L2*theta'))+(Kr*Xr) ;

% Vr'=(1/Mr)[-Ksr(Xr+(L2*theta))-Bsr(Vr+(L2*theta'))+(Kr*Xr)] ;

% MbXb"=Ksf(Xf-[Xb-(L1*theta)])+Bsf(Xf'-[Xb'-(L1*theta')])+Ksr(Xr+[Xb+(L2*theta)])+Bsr(Xr'+[Xb'+(L2*theta'))+fa(t);

% Vb'=(1/Mb)[Ksf(Xf-(L1*theta))+Bsf(Vf-(L1*omega))+Ksr(Xr+(L2*theta))+Bsr(Vr+(L2*omega))+fa(t)];

syms Ksf Xf L1 theta Bsf Vf omega Ksr Bsr Kf Kr Vr Mr Ar Xr L2 Af Ab Mf Mb fa(t) Alpha Ic L3 Xb Vb

z=Mr*Ar==-Ksr*(Xr+(Xb+(L2*theta)))-Bsr*(Vr+(Vb+(L2*omega)))+(Kr*Xr);

expand(solve(z,Ar))

ans =

y=Mf*Af==-Ksf*(Xf-(Xb-(L1*theta)))-Bsf*(Vf-(Vb-(L1*omega)))+(Kf*Xf);

expand(solve(y,Af))

ans =

x=Mb*Ab==Ksf*(Xf-(Xb-(L1*theta)))+Bsf*(Vf-(Vb-(L1*omega)))+Ksr*(Xr+(Xb+(L2*theta)))+Bsr*(Vr+(Vb+(L2*omega)))+fa(t);

expand(solve(x,Ab))

ans =

%.....................................................................SUM_OF_EXTERNAL_TORQUE=J*THETA"...............................................................................................

% Ic*theta"={-[Ksf(Xf-(L1*theta))*L1]-[Bsf(Xf'-(L1*theta'))*L1]+[Ksr(Xr+(L2*theta))*L2]+[Bsr(Xr'+(L2*theta'))*L2]+[fa(t)*L3]};

% omega'=(1/Ic){[-Ksf*Xf + Ksf((L1)^2)*theta)] + [-Bsf*Vf*L1 + Bsf*Vf*L1 + Bsf((L1)^2)*omega)] + [Ksr*Xr*L2 + Ksr((L2)^2)*theta)] + [Bsr*Vr*L2 + Bsr((L2)^2)*omega)*L2] + [fa(t)*L3]};

w=Ic*Alpha==-Ksf*(Xf-(Xb-(L1*theta)))*L1-Bsf*(Vf-(Vb-(L1*omega)))*L1+Ksr*(Xr+(Xb+(L2*theta)))*L2+Bsr*(Vr+(Vb+(L2*omega)))*L2+fa(t)*L3;

expand(solve(w,Alpha))

ans =

%.....................................................................OUTPUT_EQUATIONS..................................................................................................................

syms f_Ksf f_Kf f_Bsf f_Ksr f_Kr f_Bsr

% f_Ksf=Ksf(Xf-(-L1*theta));

% f_Kf=Kf*Xf;

% f_Bsf=Bsf(Vf-(-L1*omega));

% f_Ksr=Ksr(Xr-(L2*theta));

% f_Kr=Kr*Xr;

% f_Bsr=Bsr(Vr-(L2*omega));

f_Ksf=Ksf*(Xf-(Xb-(-L1*theta)));

A=simplify(f_Ksf);

expand(A)

ans =

f_Bsf=Bsf*(Vf-(Vb-(-L1*omega)));

B=simplify(f_Bsf);

expand(B)

ans =

f_Ksr=Ksr*(Xr-(Xb+(L2*theta)));

C=simplify(f_Ksr);

expand(C)

ans =

f_Bsr=Bsr*(Vr-(Vb+(L2*omega)));

D=simplify(f_Bsr);

expand(D)

ans =

% Mb*Vb'-Ksf(Xf-(L1*theta))-Bsf(Xf'-(L1*theta'))-Ksr(Xr+(L2*theta))-Bsr(Xr'+(L2*theta'))=fa(t);

% Ic*omega'+[Ksf(Xf-(L1*theta))*L1]+[Bsf(Vf-(L1*omega))*L1]-[Ksr(Xr+(L2*theta))*L2]-[Bsr(Vr+(L2*omega))*L2]=[fa(t)*L3]

clc

clear

%-------------------------------------------------SYSTEM_PARAMETERS----------------------------------------------------------------------------------------------------------------------------------------------------------

Ic=1356; %kg-m^2

Mb=730; %kg

Mf=59; %kg

Mr=45; %kg

Kf=23000; %N/m

Ksf=18750; %N/m

Kr=16182; %N/m

Ksr=12574; %N/m

Bsf=100; %N*s/m

Bsr=100; %N*s/m

L1=1.45; %m

L2=1.39; %m

L3=0.67; %m

%-------------------------------------------------------------------------------------Matrix's-----------------------------------------------------------------------------------------------------------------------------------

%__________________________________________________INPUT__MATRIX_________________________________________________________________________________________________________________________________________________________________

%________________________________________________________________________________________________________________________________________________________________________________________________________________________________

%......Xf...................Vf......................Xr...................Vr..................Xb...........................Vb................................theta.............................omega..............................

%________________________________________________________________________________________________________________________________________________________________________________________________________________________________

a=[ 0 1 0 0 0 0 0 0; %Xf'

((-Ksf+Kf)/Mf) (Bsf/Mf) 0 0 (Ksf/Mf) (Bsf/Mf) ((Ksf*L1)/Mf) ((Bsf*L1)/Mf); %Vf'

0 0 0 1 0 0 0 0; %Xr'

0 0 ((-Ksr+Kr)/Mr) (-Bsr/Mr) (-Ksr/Mr) (-Bsr/Mr) ((-Ksr*L2)/Mr) ((-Bsr*L2)/Mr); %Vr'

0 0 0 0 0 1 0 0; %Xb'

(Ksf/Mb) (Bsf/Mb) (Ksr/Mb) (Bsr/Mb) ((Ksr-Ksf)/Mb) ((Bsr-Bsf)/Mb) (((-Ksf*L1)+(Ksr*L2))/Mb) (((-Bsf*L1)+(Bsr*L2))/Mb); %Vb'

0 0 0 0 0 0 0 1; %theta'

((-Ksf*L1)/Ic) ((-Bsf*L1)/Ic) ((Ksr*L2)/Ic) ((Bsr*L2)/Ic) (((Ksf*L1)+(Ksr*L2))/Ic) (((Bsf*L1)+(Bsr*L2))/Ic) (((Ksr*(L2^2))+(Ksf*(L1^2)))/Ic) (((Bsr*(L2^2))+(Bsf*(L1^2)))/Ic)]; %omega'

b=[0;

0;

0;

0;

0;

(1/Mb);

0;

((L3)/Ic);];

%________________________________________________OUTPUT__MATRIX_____________________________________________________________________________________________________________________________________________________________________

c=[ (Ksf) 0 0 0 (-Ksf) 0 (-Ksf*L1) 0; %f_Ksf

Kf 0 0 0 0 0 0 0; %f_Kf

0 (Bsf) 0 0 0 (-Bsf) 0 (Bsf*L1); %f_Bsf

0 0 (Ksr) 0 (-Ksr) 0 (-Ksr*L2) 0; %f_Ksr

0 0 Kr 0 0 0 0 0; %f_Kr

0 0 0 (Bsr) 0 (-Bsr) 0 (-Bsr*L2); %f_Bsr

(Ksf/Mb) (Bsf/Mb) (Ksr/Mb) (Bsr/Mb) 0 0 (((-Ksf*L1)+(Ksr*L2))/Mb) (((-Bsf*L1)+(Bsr*L2))/Mb); %Ab

0 0 0 0 (((Ksr*L2)-(Ksf*L1))/Ic) (((Bsr*L2)-(Bsf*L1))/Ic) (((Ksr*(L2^2))-(Ksf*(L1^2)))/Ic) (((Bsr*(L2^2))-(Bsf*(L1^2)))/Ic)]; %ALPHAm

d=[0;

0;

0;

0;

0;

0;

1/Mb;

(L3/Ic)];

CSTR=ss(a,b,c,d)

CSTR =

A =

x1 x2 x3 x4 x5 x6 x7 x8

x1 0 1 0 0 0 0 0 0

x2 72.03 1.695 0 0 317.8 1.695 460.8 2.458

x3 0 0 0 1 0 0 0 0

x4 0 0 80.18 -2.222 -279.4 -2.222 -388.4 -3.089

x5 0 0 0 0 0 1 0 0

x6 25.68 0.137 17.22 0.137 -8.46 0 -13.3 -0.008219

x7 0 0 0 0 0 0 0 1

x8 -20.05 -0.1069 12.89 0.1025 32.94 0.2094 46.99 0.2975

B =

u1

x1 0

x2 0

x3 0

x4 0

x5 0

x6 0.00137

x7 0

x8 0.0004941

C =

x1 x2 x3 x4 x5 x6 x7 x8

y1 1.875e+04 0 0 0 -1.875e+04 0 -2.719e+04 0

y2 2.3e+04 0 0 0 0 0 0 0

y3 0 100 0 0 0 -100 0 145

y4 0 0 1.257e+04 0 -1.257e+04 0 -1.748e+04 0

y5 0 0 1.618e+04 0 0 0 0 0

y6 0 0 0 100 0 -100 0 -139

y7 25.68 0.137 17.22 0.137 0 0 -13.3 -0.008219

y8 0 0 0 0 -7.161 -0.004425 -11.16 -0.01257

D =

u1

y1 0

y2 0

y3 0

y4 0

y5 0

y6 0

y7 0.00137

y8 0.0004941

Continuous-time state-space model.

Model Properties

sys_tf=tf(CSTR);%System transfer function

impulse(sys_tf)

sys_zpg=zpk(CSTR);%system zero pole gain

impulse(sys_zpg)

syms s real;

i=eye(8);

h=simplify(c*inv((s*i)-a)*b+d)

h =

h(1)

ans =

simplify(h(1))

ans =

%-------------------------SIMULATION-----------------------------------------------

sim('Pagels_Jon_Project_Matrix_simulation')

fat=10*(exp(-5*t));

Xf=q(:,1);

Vf=q(:,2);

Xr=q(:,3);

Vr=q(:,4);

Xb=q(:,5);

Vb=q(:,6);

theta=q(:,7);

omega=q(:,8);

%----------------------------FORCES------------------------------------------------

KSF=y(:,1);

KF=y(:,2);

BSF=y(:,3);

KSR=y(:,4);

KR=y(:,5);

BSR=y(:,6);

Ab=y(:,7);

Alpha_m=y(:,8);

%----------------------------FIGURES-------------------------------------------------

figure(1)

set(gcf,'color','white')

subplot(3,1,1);

plot(t,Xb,'-b','linewidth',0.5);grid

title('Displacement of Body Linear vs. Time')

xlabel ('Time (s)') ;

ylabel ('Displacement (m)');

subplot(3,1,2);

plot(t,theta,'-k','linewidth',0.5);grid

title('Displacement of Body Angular vs. Time ')

xlabel ('Time (s)') ;

ylabel ('Displacement (rad)');

subplot(3,1,3);

plot(t,fat,'-m','linewidth',0.5);grid

title('Impulse Force from Road vs. Time')

xlabel ('Time (s)') ;

ylabel ('Force (N/s)');

figure(2)

set(gcf,'color','white')

subplot(3,1,1);

plot(t,Xf,'-c','linewidth',0.5);grid

title('Displacement of Front Axel vs. Time')

xlabel ('Time (s)') ;

ylabel ('Displacement (m)');

subplot(3,1,2);

plot(t,KSF,'-m','linewidth',0.5);grid

hold on

subplot(3,1,2);

plot(t,KF,'-y','linewidth',0.5);grid

title('Front Spring Force vs. Time ')

xlabel ('Time (s)') ;

ylabel ('Spring Force-Front (N)');

hold off

subplot(3,1,3);

plot(t,BSF,'-g','linewidth',0.5);grid

title('Front Damper Force vs. Time')

xlabel ('Time (s)') ;

ylabel ('Damper Force-Front (N/s)');

figure(3)

set(gcf,'color','white')

subplot(3,1,1);

plot(t,Xr,'-c','linewidth',0.5);grid

title('Displacement of Rear Axel vs. Time')

xlabel ('Time (s)') ;

ylabel ('Displacement (m)');

subplot(3,1,2);

plot(t,KSR,'-m','linewidth',0.5);grid

hold on

subplot(3,1,2);

plot(t,KR,'-y','linewidth',0.5);grid

title('Rear Spring Force vs. Time ')

xlabel ('Time (s)') ;

ylabel ('Spring Force-Rear (N)');

hold off

subplot(3,1,3);

plot(t,BSR,'-g','linewidth',0.5);grid

title('Rear Damper Force vs. Time')

xlabel ('Time (s)') ;

ylabel ('Damper Force-Rear (N/s)');

I am not seeing the response in the springs and dampers everything plots as an exponential curve, I'd expect to see some oscilation. I'm not sure if I'm entering something incorrectly.