Find Q and R matriks to get eigenvalue positive

My code

clear; close; clc;

syms a1_head a2_head b hstar 

%Parameter Massa
m1 = 8095; % massa train set 1 dalam kg
m2 = 8500; % massa train set 2 dalam kg
g  = 10;

c_0_1 = 0.01176;
c_1_1 = 0.00077616;
c_2_1 = 4.48 ;

c_0_2 = 0.01176 ;
c_1_2 = 0.00077616;
c_2_2 = 4.48;

v_0 = 300;

hstar = 120;
a2_head
a_1   = -1./m1.*(c_1_1 + 2.*c_2_1.*v_0);  
a_2   = -1./m2.*(c_1_2 + 2.*c_2_2.*v_0); 


a_1_head = 1-(a_1.*hstar);
a_2_head = 1-(a_2.*hstar);
b = 1;

% Model data
A = sym(zeros(4,4));
A(1,2) = a_1_head;
A(3,2) = (a_2_head) - 1;  A(3,4) = a_2_head;
display(A);

B = sym(zeros(4,2));
B(1,1) = -b*hstar;
B(2,1) = b;
B(3,2) = -b*hstar   ;
B(4,1) = -b; B(4,2) = b;
display(B);

% Q and R matrices for ARE
Q = sym(zeros(4,4)); Q(1,:) = [10000 100 100 100]; Q(2,:) = [100 10000 100 100]; Q(3,:) = [100 100 10000 100]; Q(4,:) = [100 100 100 10000]; display(Q);
R = sym(zeros(2,2)); R(1,:) = [2500 90]; R(2,:) = [90 5000]; display(R);

% S matrix Value
s1_1 =  -(124919^(1/2)*50000000^(1/2))/60000;
s1_2 =   (124919^(1/2)*50000000^(1/2))/60000;

s2_1 = -5.2761e+03;
s2_2 =  5.2761e+03;

s3_1 = 13.4620;
s3_2 = 16.1414;

s4_1 = -232.2135;
s4_2 = 322.5401;

s7_1 = -1.5834e+05;
s7_2 =  1.5522e+05;

s12_1 = 19.1007;
s12_2 = 4.8186;

S_matrix_1 = [s1_1 s2_1 s3_1  s4_1;
              s2_1 s1_1 s7_1  s2_1;
              s3_1 s7_1 s1_1 s12_1;
              s4_1 s2_1 s12_1 s1_1
            ];
          
S_matrix_2 = [s1_2 s2_2 s3_2  s4_2;
              s2_2 s1_2 s7_2  s2_2;
              s3_2 s7_2 s1_2 s12_2;
              s4_2 s2_2 s12_2 s1_2
            ];
        
eig_1 = eig(S_matrix_1);
eig_2 = eig(S_matrix_2);

K1 = inv(R)*transpose(B)*S_matrix_1;
display(K1);
K2 = inv(R)*transpose(B)*S_matrix_2;
display(K2);

% Matrix S to find
svar = sym('s',[1 16]);
S = [svar(1:4); svar(5:8); svar(9:12); svar(13:16)];
S(2,1) = svar(2);
S(2,2) = svar(1);
S(2,4) = svar(2);
S(3,1) = svar(3);
S(3,2) = svar(7);
S(3,3) = svar(1);
S(4,1) = svar(4);
S(4,2) = svar(2);
S(4,3) = svar(12);
S(4,4) = svar(1);
display(S);

% LHS of ARE: A'*S + S*A' - S*B*Rinv*B'*S
left_ARE = transpose(A)*S + S*A - S*B*inv(R)*transpose(B)*S; 
display(left_ARE);

% RHS of ARE: A'*S + S*A' - S*B*Rinv*B'*S
right_ARE = -Q;
display(right_ARE);

% Find S in ARE
C = left_ARE;
display(C);

D = right_ARE;
display(D);

% Element Matriks S
C1 = C(1,1);
display(C1);

C2 = C(1,2);
display(C2);

C3 = C(1,3);
display(C3);

C4 = C(1,4);
display(C4);

C5 = C(2,1);
display(C5);

C6 = C(2,2);
display(C6);

C7 = C(2,3);
display(C7);

C8 = C(2,4);
display(C8);

C9 = C(3,1);
display(C9);

C10 = C(3,2);
display(C10);

C11 = C(3,3);
display(C11);

C12 = C(3,4);
display(C12);

C13 = C(4,1);
display(C13);

C14 = C(4,2);
display(C14);

C15 = C(4,3);
display(C15);

C16 = C(4,4);
display(C16);

D1 = D(1,1);
display(D1);

D2 = D(1,2);
display(D2);

D3 = D(1,3);
display(D3);

D4 = D(1,4);
display(D4);

D5 = D(2,1);
display(D5);

D6 = D(2,2);
display(D6);

D7 = D(2,3);
display(D7);

D8 = D(2,4);
display(D8);

D9 = D(3,1);
display(D9);

D10 = D(3,2);
display(D10);

D11 = D(3,3);
display(D11);

D12 = D(3,4);
display(D12);

D13 = D(4,1);
display(D13);

D14 = D(4,2);
display(D14);

D15 = D(4,3);
display(D15);

D16 = D(4,4);
display(D16);

E1 = C1 == D1;
display(E1);

E2 = C2 == D2;
display(E2);

E3 = C3 == D3;
display(E3)

E4 = C4 == D4;
display(E4);

E5 = C5 == D5;
display(E5);

E6 = C6 == D6;
display(E6);

E7 = C7 == D7;
display(E7);

E8 = C8 == D8;
display(E8);

E9 = C9 == D9;
display(E9);

E10 = C10 == D10;
display(E10);

E11 = C11 == D11;
display(E11);

E12 = C12 == D12;
display(E12);

E13 = C13 == D13;
display(E13);

E14 = C14 == D14;
display(E14);

E15 = C15 == D15;
display(E15);

E16 = C16 == D16;
display(E16);

syms s1 s2 s3 s4 s5 s6 s7 s8 s9 s10 s11 s12 s13 s14 s15 s16
[Sol_s1_1, Sol_s2_1, Sol_s3_1, Sol_s4_1] = solve(E1,s1,s2,s3,s4)

[Sol_s2_2, Sol_s3_2, Sol_s4_2, Sol_s7_1] = solve(E2,s2,s3,s4,s7,'Real',true)

[Sol_s2_3, Sol_s3_3, Sol_s4_3, Sol_s7_2, sol_s12_1] = solve(E3,s2,s3,s4,s7,s12,'Real',true) 

[Sol_s3_4, Sol_s4_4,sol_s12_2] = solve(E4,s3,s4,s12,'Real',true) 

[Sol_s4_5, Sol_s7_3] = solve(E5,s4,s7,'Real',true)

[Sol_s7_4] = solve(E6,s7,'Real',true)
  
[Sol_s7_5, Sol_s12_3] = solve(E7,s7,s12,'Real',true)

[Sol_s12_4] = solve(E8,s12,'Real',true)
% 
% [Sol_s2_9, Sol_s3_7, Sol_s4_7, Sol_s7_7, Sol_s12_5] = solve(E9,s2,s3,s4,s7,s12,'Real',true)
% 
% [Sol_s3_8, Sol_s7_8, Sol_s12_6] = solve(E10,s3,s7,s12,'Real',true)
% 
% [Sol_s3_9, Sol_s7_9, Sol_s12_7] = solve(E11,s3,s7,s12,'Real',true)
% 
% [Sol_s4_8, Sol_s7_10, Sol_s12_8] = solve(E12,s4,s7,s12,'Real',true)
% 
% [Sol_s4_9, Sol_s12_9] = solve(E13,s4,s12,'Real',true)
% 
% [Sol_s7_11, Sol_s12_10] = solve(E14,s7,s12,'Real',true)
% 
% [Sol_s12_11] = solve(E15,s12,'Real',true)
% 
% [Sol_s12_12] = solve(E16,s12,'Real',true)

% syms s1 s2 s3 s4 s5 s6 s7 s8 s9 s10 s11 s12 s13 s14 s15 s16
% [Sol_s1_1, Sol_s2_1, Sol_s3_1, Sol_s4_1] = solve(E1,s1,s2,s3,s4)
% 
% [Sol_s2_2, Sol_s3_2, Sol_s4_2, Sol_s7_1] = solve(E2,s2,s3,s4,s7,'Real',true)
%  
% [Sol_s3_3, Sol_s4_3, Sol_s7_2, sol_s12_1] = solve(E3,s3,s4,s7,s12,'Real',true) 
%  
% [Sol_s4_4, Sol_s12_2] = solve(E4,s4,s12,'Real',true) 
% 
% [Sol_s7_3] = solve(E5,s7,'Real',true)
% % 
% % [Sol_s8_4] = solve(E6,s8,'Real',true)
% %   
% [Sol_s12_3] = solve(E7,s12)
%  
% % [Sol_s8_6, Sol_s12_4, Sol_s16] = solve(E8,s8,s12,s16,'Real',true)
% %  
% % [Sol_s11_5, Sol_s12_5] = solve(E9,s11,s12,'Real',true)
% %  
% % [Sol_s12_6] = solve(E10,s12,'Real',true)
% %  
% % [Sol_s12_7] = solve(E11,s12,'Real',true)
% % 
% % [Sol_s16_3] = solve(E12,s16,'Real',true)
% % % % 
% % % [Sol_s8_9, Sol_s12_9, Sol_s16_4] = solve(E13,s8,s12,s16,'Real',true)
% % % 
% % % [Sol_s2_22, Sol_s4_10, Sol_s6_7, Sol_s7_11, Sol_s8_10, Sol_s12_10, Sol_s16_5] = solve(E14,s2,s4,s6,s7,s8,s12,s16,'Real',true)
% % % 
% % % [Sol_s3_12, Sol_s4_11, Sol_s7_12, Sol_s8_11, Sol_s11_9, Sol_s12_11, Sol_s16_6] = solve(E15,s3,s4,s7,s8,s11,s12,s16,'Real',true)
% % % 
% % % [Sol_s4_12, Sol_s8_12, Sol_s12_12, Sol_s16_7] = solve(E16,s4,s8,s12,s16,'Real',true)

In this code, i need to find the s1,s2,s3,s4,s7,s12 to get my matrix S.

After i get the S matrix, I find the eigenvalue of matrix 4x4

I have problem to get the true Q and R matrix, if I can get the true Q and R matrix I can get the eigenvalue positive.