%COEFFSPOL_M20
% Code illustrating that solving the nonlinear system presented in Equation
% (33) of [1] reduces to evaluating a rational or polynomial equation in 
% the variable i2.
% y0=A^2*(sqrt(c8)*A^2+c7/(2*sqrt(c8))*A); %solution with positive sign in (28) of [3]
% y1=sum(c.*A.^([3:8]'));
% y2=(y1+d2*A^2+d1*A)*(y1+e0*y0+e2*A^2+e1*A)+f0*y1+g0*y0+h2*A^2+h1*A;
% y3=y2*(y0+i2*A^2+i1*A)+j0*y0+j2*A^2+j1*A+k0I;
%   
% Here, y3 is a polynomial of degree m = 20 in A, denoted as:
% P_m(A) = sum_{i=0}^{m} b_i * A^i, m = 20, 
% where b_i are given real or complex scalar coefficients
% Hence, y3 = P_m(A)
%
% The total cost of the above polynomial evaluation is 5 matrix multiplications.
%
% Note: This code is for illustrative purposes only. For numerical solutions,
% it is significantly more efficient to use `vpasolve`.
%
%   References: 
%
%   [1] J. Sastre, J. Ibanez, J.M. alonso E. Defez, Beyond 
%       Paterson–Stockmeyer: Advancing Matrix Polynomial Computation, 
%       presented at the 5th Int. Conf. on Applied Mathematics, Computational
%       Science and Systems Engineering, Paris, France, Apr. 14–16, 2025, 
%       submitted to WSEAS Transactions on Mathematics.
%   [2] J. Sastre, J. Ibanez, Efficient Evaluation of Matrix Polynomials
%       beyond the Paterson–Stockmeyer Method. Mathematics 2021, 9, 1600.
%       https://doi.org/10.3390/math9141600.
%   [3] J. Sastre, Efficient evaluation of matrix polynomials, Linear 
%       Algebra Appl., 539, 2018, 229-250.
%   [4] J. Sastre, J. Ibanez, E. Defez, Boosting the computation of
%       the matrix exponential, Appl. Math. Comput. 340 (2019), 206-220.

%   Author: Jorge Sastre 2025/4/10
%   Revised version: 2025/4/10.
%
%   Group of High Performance Scientific Computing (HiPerSC)
%   Universitat Politecnica de Valencia (Spain)
%   http://hipersc.blogs.upv.es
%   Work supported by Generalitat Valenciana Grant CIAICO/2023/275

diary('coeffspol_m20.txt')

echo on

syms A c3 c4 c5 c6 c7 c8 d0 d1 d2 e0 e1 e2 f0 g0 h2 h1 i2 i1 i0 j0 j2 j1 k0I

syms b0 b1 b2 b3 b4 b5 b6 b7 b8 b9 b10 b11 b12 b13 b14 b15 b16 b17 b18 b19 b20

disp('The polynomial to evaluate is')
disp('P_m(A) = sum_{i=0}^{m} b_i * A^i, m = 20')

m = 20

c=[c3; c4; c5; c6; c7; c8];

% Polynomial P_20(A) to evaluate, defined by its coefficients (real or complex)
b=[b20; b19; b18; b17; b16; b15; b14; b13; b12; b11; b10; b9; b8; b7; b6; b5; b4; b3; b2; b1; b0];

%% Evaluation formulas
disp('Evaluation formulas more efficient than the Paterson-Stockmeyer method')
y0=A^2*(sqrt(c8)*A^2+c7/(2*sqrt(c8))*A); %solution with positive sign in (28) of [3]

y1=sum(c.*A.^([3:8]'));

y2=(y1+d2*A^2+d1*A)*(y1+e0*y0+e2*A^2+e1*A)+f0*y1+g0*y0+h2*A^2+h1*A;

y3=y2*(y0+i2*A^2+i1*A)+j0*y0+j2*A^2+j1*A+k0I;

%% Equations
disp('System of m+1 = 21 equations to solve')
[cy3,a1]=coeffs(y3,A);

cy3=cy3.';

v=[cy3 b a1.'] %Shows the coefficients of each power of A
% [                                                                                                                                                                                                                                                                                                                                                                                                                                                 c8^(5/2), b20, A^20]
% [                                                                                                                                                                                                                                                                                                                                                                                                                                        (5*c7*c8^(3/2))/2, b19, A^19]
% [                                                                                                                                                                                                                                                                                                                                                                                                      c8^(1/2)*(c7^2 + 2*c6*c8) + c8^2*i2 + c7^2*c8^(1/2), b18, A^18]
% [                                                                                                                                                                                                                                                                                                                                                                 c8^2*i1 + c8^(1/2)*(2*c5*c8 + 2*c6*c7) + 2*c7*c8*i2 + (c7*(c7^2 + 2*c6*c8))/(2*c8^(1/2)), b17, A^17]
% [                                                                                                                                                                                                                                                                                                                     i2*(c7^2 + 2*c6*c8) + c8^(1/2)*(c6^2 + c8*(c4 + c8^(1/2)*e0) + c4*c8 + 2*c5*c7) + (c7*(2*c5*c8 + 2*c6*c7))/(2*c8^(1/2)) + 2*c7*c8*i1, b16, A^16]
% [                                                                                                                                                                                                                                          i1*(c7^2 + 2*c6*c8) + c8^(1/2)*(c7*(c4 + c8^(1/2)*e0) + c3*c8 + c4*c7 + 2*c5*c6 + c8*(c3 + (c7*e0)/(2*c8^(1/2)))) + i2*(2*c5*c8 + 2*c6*c7) + (c7*(c6^2 + c8*(c4 + c8^(1/2)*e0) + c4*c8 + 2*c5*c7))/(2*c8^(1/2)), b15, A^15]
% [                                                                                                                                                           i2*(c6^2 + c8*(c4 + c8^(1/2)*e0) + c4*c8 + 2*c5*c7) + c8^(1/2)*(c5^2 + c6*(c4 + c8^(1/2)*e0) + c3*c7 + c4*c6 + c8*d2 + c8*e2 + c7*(c3 + (c7*e0)/(2*c8^(1/2)))) + i1*(2*c5*c8 + 2*c6*c7) + (c7*(c7*(c4 + c8^(1/2)*e0) + c3*c8 + c4*c7 + 2*c5*c6 + c8*(c3 + (c7*e0)/(2*c8^(1/2)))))/(2*c8^(1/2)), b14, A^14]
% [                                                                      i2*(c7*(c4 + c8^(1/2)*e0) + c3*c8 + c4*c7 + 2*c5*c6 + c8*(c3 + (c7*e0)/(2*c8^(1/2)))) + i1*(c6^2 + c8*(c4 + c8^(1/2)*e0) + c4*c8 + 2*c5*c7) + c8^(1/2)*(c5*(c4 + c8^(1/2)*e0) + c3*c6 + c4*c5 + c7*d2 + c8*d1 + c7*e2 + c8*e1 + c6*(c3 + (c7*e0)/(2*c8^(1/2)))) + (c7*(c5^2 + c6*(c4 + c8^(1/2)*e0) + c3*c7 + c4*c6 + c8*d2 + c8*e2 + c7*(c3 + (c7*e0)/(2*c8^(1/2)))))/(2*c8^(1/2)), b13, A^13]
% [              i1*(c7*(c4 + c8^(1/2)*e0) + c3*c8 + c4*c7 + 2*c5*c6 + c8*(c3 + (c7*e0)/(2*c8^(1/2)))) + c8^(1/2)*(c4*(c4 + c8^(1/2)*e0) + c3*c5 + c6*d2 + c7*d1 + c6*e2 + c7*e1 + c8*f0 + c5*(c3 + (c7*e0)/(2*c8^(1/2)))) + i2*(c5^2 + c6*(c4 + c8^(1/2)*e0) + c3*c7 + c4*c6 + c8*d2 + c8*e2 + c7*(c3 + (c7*e0)/(2*c8^(1/2)))) + (c7*(c5*(c4 + c8^(1/2)*e0) + c3*c6 + c4*c5 + c7*d2 + c8*d1 + c7*e2 + c8*e1 + c6*(c3 + (c7*e0)/(2*c8^(1/2)))))/(2*c8^(1/2)), b12, A^12]
% [i2*(c5*(c4 + c8^(1/2)*e0) + c3*c6 + c4*c5 + c7*d2 + c8*d1 + c7*e2 + c8*e1 + c6*(c3 + (c7*e0)/(2*c8^(1/2)))) + c8^(1/2)*(c3*(c4 + c8^(1/2)*e0) + c5*d2 + c6*d1 + c5*e2 + c6*e1 + c7*f0 + c4*(c3 + (c7*e0)/(2*c8^(1/2)))) + i1*(c5^2 + c6*(c4 + c8^(1/2)*e0) + c3*c7 + c4*c6 + c8*d2 + c8*e2 + c7*(c3 + (c7*e0)/(2*c8^(1/2)))) + (c7*(c4*(c4 + c8^(1/2)*e0) + c3*c5 + c6*d2 + c7*d1 + c6*e2 + c7*e1 + c8*f0 + c5*(c3 + (c7*e0)/(2*c8^(1/2)))))/(2*c8^(1/2)), b11, A^11]
% [       i1*(c5*(c4 + c8^(1/2)*e0) + c3*c6 + c4*c5 + c7*d2 + c8*d1 + c7*e2 + c8*e1 + c6*(c3 + (c7*e0)/(2*c8^(1/2)))) + i2*(c4*(c4 + c8^(1/2)*e0) + c3*c5 + c6*d2 + c7*d1 + c6*e2 + c7*e1 + c8*f0 + c5*(c3 + (c7*e0)/(2*c8^(1/2)))) + c8^(1/2)*(d2*(c4 + c8^(1/2)*e0) + c5*d1 + c4*e2 + c5*e1 + c6*f0 + c3*(c3 + (c7*e0)/(2*c8^(1/2)))) + (c7*(c3*(c4 + c8^(1/2)*e0) + c5*d2 + c6*d1 + c5*e2 + c6*e1 + c7*f0 + c4*(c3 + (c7*e0)/(2*c8^(1/2)))))/(2*c8^(1/2)), b10, A^10]
% [                               i2*(c3*(c4 + c8^(1/2)*e0) + c5*d2 + c6*d1 + c5*e2 + c6*e1 + c7*f0 + c4*(c3 + (c7*e0)/(2*c8^(1/2)))) + c8^(1/2)*(d1*(c4 + c8^(1/2)*e0) + c3*e2 + c4*e1 + c5*f0 + d2*(c3 + (c7*e0)/(2*c8^(1/2)))) + i1*(c4*(c4 + c8^(1/2)*e0) + c3*c5 + c6*d2 + c7*d1 + c6*e2 + c7*e1 + c8*f0 + c5*(c3 + (c7*e0)/(2*c8^(1/2)))) + (c7*(d2*(c4 + c8^(1/2)*e0) + c5*d1 + c4*e2 + c5*e1 + c6*f0 + c3*(c3 + (c7*e0)/(2*c8^(1/2)))))/(2*c8^(1/2)),  b9,  A^9]
% [                                                                 i2*(d2*(c4 + c8^(1/2)*e0) + c5*d1 + c4*e2 + c5*e1 + c6*f0 + c3*(c3 + (c7*e0)/(2*c8^(1/2)))) + i1*(c3*(c4 + c8^(1/2)*e0) + c5*d2 + c6*d1 + c5*e2 + c6*e1 + c7*f0 + c4*(c3 + (c7*e0)/(2*c8^(1/2)))) + c8^(1/2)*(c3*e1 + c4*f0 + d2*e2 + d1*(c3 + (c7*e0)/(2*c8^(1/2))) + c8^(1/2)*g0) + (c7*(d1*(c4 + c8^(1/2)*e0) + c3*e2 + c4*e1 + c5*f0 + d2*(c3 + (c7*e0)/(2*c8^(1/2)))))/(2*c8^(1/2)),  b8,  A^8]
% [                                                                                                                   i2*(d1*(c4 + c8^(1/2)*e0) + c3*e2 + c4*e1 + c5*f0 + d2*(c3 + (c7*e0)/(2*c8^(1/2)))) + i1*(d2*(c4 + c8^(1/2)*e0) + c5*d1 + c4*e2 + c5*e1 + c6*f0 + c3*(c3 + (c7*e0)/(2*c8^(1/2)))) + c8^(1/2)*(c3*f0 + d1*e2 + d2*e1 + (c7*g0)/(2*c8^(1/2))) + (c7*(c3*e1 + c4*f0 + d2*e2 + d1*(c3 + (c7*e0)/(2*c8^(1/2))) + c8^(1/2)*g0))/(2*c8^(1/2)),  b7,  A^7]
% [                                                                                                                                                                                               c8^(1/2)*(h2 + d1*e1) + i1*(d1*(c4 + c8^(1/2)*e0) + c3*e2 + c4*e1 + c5*f0 + d2*(c3 + (c7*e0)/(2*c8^(1/2)))) + i2*(c3*e1 + c4*f0 + d2*e2 + d1*(c3 + (c7*e0)/(2*c8^(1/2))) + c8^(1/2)*g0) + (c7*(c3*f0 + d1*e2 + d2*e1 + (c7*g0)/(2*c8^(1/2))))/(2*c8^(1/2)),  b6,  A^6]
% [                                                                                                                                                                                                                                                                             c8^(1/2)*h1 + i2*(c3*f0 + d1*e2 + d2*e1 + (c7*g0)/(2*c8^(1/2))) + i1*(c3*e1 + c4*f0 + d2*e2 + d1*(c3 + (c7*e0)/(2*c8^(1/2))) + c8^(1/2)*g0) + (c7*(h2 + d1*e1))/(2*c8^(1/2)),  b5,  A^5]
% [                                                                                                                                                                                                                                                                                                                                                 c8^(1/2)*j0 + i1*(c3*f0 + d1*e2 + d2*e1 + (c7*g0)/(2*c8^(1/2))) + i2*(h2 + d1*e1) + (c7*h1)/(2*c8^(1/2)),  b4,  A^4]
% [                                                                                                                                                                                                                                                                                                                                                                                                           h1*i2 + i1*(h2 + d1*e1) + (c7*j0)/(2*c8^(1/2)),  b3,  A^3]
% [                                                                                                                                                                                                                                                                                                                                                                                                                                               j2 + h1*i1,  b2,  A^2]
% [                                                                                                                                                                                                                                                                                                                                                                                                                                                       j1,  b1,    A]
% [                                                                                                                                                                                                                                                                                                                                                                                                                                                      k0I,  b0,    1]

%% System of 21 equations

cy3=cy3-b; 

%% Variable solving and substitution

disp('Variable solving and substitution')

c8s=solve(cy3(1),c8) %3 solutions

c8s=c8s(1) %real solution

cy3=subs(cy3,c8,c8s); %variable substitution

c7s=solve(cy3(2),c7,'ReturnConditions',true); %one solution depending on c8

c7s.conditions %conditions for the existence of solutions: b20 ~= 0

c7s=c7s.c7 %solution

cy3=subs(cy3,c7,c7s); %variable substitution

symvar(cy3(3))  %depends on the unknown variables [c6, i2] and the polynomial coefficients [b18, b19, b20]

c6s=solve(cy3(3), c6,'ReturnConditions',true);

c6s.conditions %Conditions for the existence of solutions: b20 ~= 0

c6s=c6s.c6 %one solution depending on i2 and the polynomial coefficients [b18, b19, b20]

cy3=subs(cy3,c6,c6s); %variable substitution

symvar(cy3(4)) %depends on the unknown variables [c5, i1, i2]  and the polynomial coefficients [b17, b18, b19, b20]

c5s=solve(cy3(4), c5,'ReturnConditions',true); %one solution depending on the unknown variables [i1, i2] and the polynomial coefficients [b17, b18, b19, b20]
%since i1 depends on i2 (see below), c5s depends on %one solution depending on i2 and the polynomial coefficients [b17, b18, b19, b20]

c5s.conditions %conditions for the existence of solutions: b20 ~= 0

c5s=c5s.c5 %one solution depending on i2 and the polynomial coefficients [b17, b18, b19, b20]

cy3=simplify(subs(cy3,c5,c5s)); %variable substitution

symvar(cy3(5))  %depends on [c4, e0, i1, i2] and [b16, b17, b18, b19, b20]

e0s=solve(cy3(5), e0,'ReturnConditions',true); %depends on [c4, i1, i2] and [b16, b17, b18, b19, b20]
%since c4 and i1 depend on i2 (see below) then e0s depends on i2 and [b16, b17, b18, b19, b20]

e0s.conditions %conditions: b20 ~= 0

e0s=e0s.e0; %one solution depending on i2 and [b16, b17, b18, b19, b20]

cy3=simplify(subs(cy3,e0,e0s)); %variable substitution

symvar(cy3(6))  %depends on [c3, c4, i1, i2] and [b15, b16, b17, b18, b19, b20]

c3s=solve(cy3(6), c3,'ReturnConditions',true); %depends on [c4, i1, i2] and [b15, b16, b17, b18, b19, b20]
%since c4 and i1 depend on i2 (see below) then c3s depends on i2 and [b15, b16, b17, b18, b19, b20]

c3s.conditions %conditions: b20 ~= 0

c3s=c3s.c3; %one solution depending on i2 and [b15, b16, b17, b18, b19, b20]

cy3=simplify(subs(cy3,c3,c3s));

symvar(cy3(7)) %depends on [d2, e2, i1, i2] and [b14, b15, b16, b17, b18, b19, b20]

d2s=solve(cy3(7), d2,'ReturnConditions',true); %depends on [e2, i1, i2] and [b14, b15, b16, b17, b18, b19, b20]
%since e2 and i1 depend on i2 (see below) then d2s depends on i2 and [b14, b15, b16, b17, b18, b19, b20]

d2s.conditions %b20 ~= 0

d2s=d2s.d2; %one solution depending on i2 and [b14, b15, b16, b17, b18, b19, b20]

cy3=simplify(subs(cy3,d2,d2s));

symvar(cy3(8)) %[d1, e1, i1, i2] and [b13, b14, b15, b16, b17, b18, b19, b20]

d1s=solve(cy3(8), d1,'ReturnConditions',true); %depends on %[e1, i1, i2] and [b13, b14, b15, b16, b17, b18, b19, b20]
%since e1 and i1 depend on i2 (see below) then d1s depends on i2 and [b13, b14, b15, b16, b17, b18, b19, b20]

d1s.conditions %b20 ~= 0

d1s=d1s.d1; %one solution depending on i2 and [b13, b14, b15, b16, b17, b18, b19, b20]

cy3=simplify(subs(cy3,d1,d1s));

symvar(cy3(9)) %[c4, f0, i1, i2] and [b12, b13, b14, b15, b16, b17, b18, b19, b20]

f0s=solve(cy3(9), f0,'ReturnConditions',true); %depends on [c4, i1, i2] and [b12, b13, b14, b15, b16, b17, b18, b19, b20]
%since i1 and c4 depend on i2 (see below) then f0s depends on i2 and 
% [b12, b13, b14, b15, b16, b17, b18, b19, b20]

f0s.conditions %b20 ~= 0

f0s=f0s.f0; %one solution depending on i2 and [b12, b13, b14, b15, b16, b17, b18, b19, b20]

cy3=simplify(subs(cy3,f0,f0s));

symvar(cy3(10)) %Depends on [i1, i2] and [b11, b12, b13, b14, b15, b16, b17, b18, b19, b20]
[cc,aa]=coeffs(cy3(10),i1) %[i1^3, i1^2, i1, 1] polynomial of degree 3 on variable i1 with coefficients depending on i2 and the polynomial coefficients to evaluate
[cc,aa]=coeffs(cy3(10),i2) %[i2^4, i2^3, i2^2, i2, 1] polynomial of degree 4 on variable i2 with coefficients depending on i1 and the polynomial coefficients to evaluate
i1s=solve(cy3(10), i1,'ReturnConditions',true)
%i1 is the solution of a cubic polynomial equation whose coefficients 
% depend on i2 and the coefficients bi of the polynomail to evaluate, then 
% there are three solutions of i1 depending on i2 and the polynomial coefficients

i1s.conditions %b20 ~= 0

symvar(cy3(11)) %depends on [c4, e2, i1, i2] and [b10, b12, b13, b14, b15, b16, b17, b18, b19, b20]

e2s=solve(cy3(11), e2,'ReturnConditions',true); %depends on [c4, i1, i2] and [b10, b12, b13, b14, b15, b16, b17, b18, b19, b20]
%since i1 depends on i2 and c4 also depends on i2 (see below) then e2s
%depends on i2 and [b10, b12, b13, b14, b15, b16, b17, b18, b19, b20]

e2s.conditions %500*b16*b20^(43/5) + 375*b20^(46/5)*i2^2 + 220*b18*b19^2*b20^(33/5) + 100*b19^2*b20^(37/5)*i2 ~= 1000*b20^(46/5)*c4 + 48*b19^4*b20^(28/5) + 125*b18^2*b20^(38/5) + 300*b17*b19*b20^(38/5) + 250*b18*b20^(42/5)*i2 + 100*b19*b20^(42/5)*i1

e2s=e2s.e2; %one solution depending on i2 and [b10, b12, b13, b14, b15, b16, b17, b18, b19, b20]

cy3=simplify(subs(cy3,e2,e2s));

symvar(cy3(12)) %depends on [c4, e1, i1, i2] and [b10, b12, b13, b14, b15, b16, b17, b18, b19, b20, b9]

e1s=solve(cy3(12), e1,'ReturnConditions',true);%depends on [c4, i1, i2] and [b10, b12, b13, b14, b15, b16, b17, b18, b19, b20, b9]
%since i1 depends on i2 and c4 also depends on i2 (see below) then e1s
%depends on i2 and [b10, b12, b13, b14, b15, b16, b17, b18, b19, b20, b9]

e1s.conditions %500*b16*b20^(48/5) + 375*b20^(51/5)*i2^2 + 220*b18*b19^2*b20^(38/5) + 100*b19^2*b20^(42/5)*i2 ~= 1000*b20^(51/5)*c4 + 48*b19^4*b20^(33/5) + 125*b18^2*b20^(43/5) + 300*b17*b19*b20^(43/5) + 250*b18*b20^(47/5)*i2 + 100*b19*b20^(47/5)*i1

e1s=e1s.e1; %one solution depending on i2 and [b10, b12, b13, b14, b15, b16, b17, b18, b19, b20, b9]

cy3=simplify(subs(cy3,e1,e1s));

symvar(cy3(13)) %depends on [c4, g0, i1, i2] and [b10, b12, b13, b14, b15, b16, b17, b18, b19, b20, b8, b9]

g0s=solve(cy3(13), g0,'ReturnConditions',true); %depends on [c4, i1, i2] and [b10, b12, b13, b14, b15, b16, b17, b18, b19, b20, b8, b9]
%since c4 and i1 depend on i2 (see below for c4) and [b10, b12, b13, b14, b15, b16, b17, b18, b19, b20, b8, b9]
%then g0s depends on i2 and [b10, b12, b13, b14, b15, b16, b17, b18, b19, b20, b8, b9]
g0s.conditions %b20 ~= 0 & 1000*b20^(18/5)*c4 + 48*b19^4 + 125*b18^2*b20^2 + 300*b17*b19*b20^2 + 250*b18*b20^(14/5)*i2 + 100*b19*b20^(14/5)*i1 ~= 500*b16*b20^3 + 375*b20^(18/5)*i2^2 + 100*b19^2*b20^(9/5)*i2 + 220*b18*b19^2*b20

g0s=g0s.g0; %one solution depending on i2 and [b10, b12, b13, b14, b15, b16, b17, b18, b19, b20, b8, b9]

cy3=simplify(subs(cy3,g0,g0s));

symvar(cy3(14)) %Depends on [c4, i1, i2, b10, b12, b13, b14, b15, b16, b17, b18, b19, b20, b7, b8, b9]

c4s=solve(cy3(14), c4, 'ReturnConditions',true); %4 roots of a polynomial of degree 4 in c4 with coefficients depending on i1, i2 and the coefficients of the polynomial to evaluate

c4s.conditions;

c4s=c4s.c4; %4 roots of a polynomial of degree 4, depending on [i1, i2] and [b10, b12, b13, b14, b15, b16, b17, b18, b19, b20, b7, b8, b9]
% since i1 depends on i2 then c4 depends on i2 and [b10, b12, b13, b14, b15, b16, b17, b18, b19, b20, b7, b8, b9]

symvar(cy3(15)) %depends on [c4, h2, i1, i2, b10, b12, b13, b14, b15, b16, b17, b18, b19, b20, b6, b8, b9]

h2s=solve(cy3(15), h2, 'ReturnConditions',true); %one solution depending on [c4, i1, i2] and [b10, b12, b13, b14, b15, b16, b17, b18, b19, b20, b6, b8, b9]
%since c4, i1, depend on i2 then, h2s depends on i2 and and [b10, b12, b13, b14, b15, b16, b17, b18, b19, b20, b6, b8, b9]

h2s.conditions; %conditions

h2s=h2s.h2; %one solution depending on i2 and the coefficients of the polynomial to evaluate [b10, b12, b13, b14, b15, b16, b17, b18, b19, b20, b6, b8, b9]

symvar(cy3(16)) %depends on [c4, h1, h2, i1, i2] amd [b10, b12, b13, b14, b15, b16, b17, b18, b19, b20, b5, b8, b9]

h1s=solve(cy3(16), h1, 'ReturnConditions',true); %one solution depending on [c4, h2, i1, i2] and [b10, b12, b13, b14, b15, b16, b17, b18, b19, b20, b5, b8, b9]
%since c4 and h2 depend on i2 and the polynomial coefficients, then h1 depends on i2 and the polynomial coefficients

h1s.conditions;

h1s=h1s.h1; %one solution depending on i2 and the coefficients of the polynomial to evaluate [b10, b12, b13, b14, b15, b16, b17, b18, b19, b20, b5, b8, b9]

symvar(cy3(17)) %depends on [c4, h1, h2, i1, i2, j0, b10, b12, b13, b14, b15, b16, b17, b18, b19, b20, b4, b8, b9]

j0s=solve(cy3(17), j0, 'ReturnConditions',true); %one solution depending on [c4, h1, h2, i1, i2] and [b10, b12, b13, b14, b15, b16, b17, b18, b19, b20, b4, b8, b9]
%since c4, h1, h2, i1 depend on i2, then j0s depends on i2 and the
%coefficients of the polynomial to evaluate
j0s.conditions;

j0s=j0s.j0; %one solution depending on i2 and the coefficients of the polynomial to evaluate

symvar(cy3(18)) %Depends on [c4, h1, h2, i1, i2, j0[ and [b10, b12, b13, b14, b15, b16, b17, b18, b19, b20, b3, b9]
% [c4, h1, h2, i1, i2, j0] depend on i2, then equation symvar(cy3(18)) can
% be written as on rational or polynomial equation in i2.

%Taking into account (cy3(19:21)):

j2 = b2 - h1s*i1s.i1(1); %h1s and i1s has three solutions i1s.i1(k), k=1,2,3 depending on i2

j1 = b1;

k0I = b0;

diary off