There should be a lot more displayed. Look at the code and you'll see that every time the gevp solver is run, it should display the iteration of the lower and upper bound.

Also, I hope you actually try to understand the code, and not just take it for granted. Put a break in the code and follow the flow in order to understand the logic. If ilmi returns two empty matrices (is that what you claim?), then what does that mean in the code? It should be clear.

When I ran the code it returned a feasible solution after roughly 10 gevp calls, where every gevp call took 10 iterations or something like that, i.e., in total something like 100 calls to sedumi, of which 10 are displayed.

It could also be a good exercise to make the gevp code more efficient by supplying inital upper bounds, after 1 iteration you know that the last ai is an upper bound

I have written this code :
***********************************
Q = eye(2);
A=[0 1;1 0];
B=[1 0]';
C=[1 15];
delta = 1e-3;

a = sdpvar (1);
F = sdpvar(1,1);
P = sdpvar(2);
%X = sdpvar(2);

X = care(A,B*B',Q);

%i = 1;
while i <= 10
RHS = [P zeros(2,1); zeros(1,2) 0];
LHS = [A'*P+P*A-X*B*B'*P-P*B*B'*X+X*B*B'*X (B'*P+F*C)' ; (B'*P+F*C) -eye(1)];

% Find optimal ai
ai = solvelocalgevp(LHS,RHS,[],a);

% Now solve with that ai to recover the last feasible solution in the
% gevp code (could have been returned from there but it is easier to
% solve again...)
% ops = sdpsettings('debug',1) ;
Constraints = [LHS <= ai*RHS,P >= 0];
solvesdp(Constraints,trace(P))

if ai <= 0
P = double(P);
F = double(F);
break
end

Pi = double(P);
if norm(X-Pi) <= delta
P = [];
F = [];
break
end
X = Pi;
end
X
Pi= double(P) % extraire la matrice solution P
norm(X-Pi)
Fi= double(F) % extraire la matrice solution F
eig(A)
eig(A+B*Fi*C)
***********************************
with solvelocalgevp
***********************************
function t_upper = solvelocalgevp(A,B,C,t)
% Résoudre min t sous A(x)<=t*B(x), B>=0, C>=0

% Trouver des limites intiaux inf et sup
%ops = sdpsettings('verbose',0); % En mettant verbose à 0, le solveur va exécuter un affichage minimal ops : variable options
ops = sdpsettings('debug',1);
sol = solvesdp([C>=0, A <= 0*B, B>=0],[],ops); % solvesdp(contraintes ,Objectif ,options)
if sol.problem == 0
% Ok, 0 est une limite sup , i.e. feasible
t_upper = 0;
t_lower = -1;
while sol.problem == 0
t_lower = (t_lower)*2;
sol = solvesdp([C>=0, A <= t_lower*B, B>=0],[],ops);
end
else sol.problem = 1;
% Ok 0 est une limite inf, i.e. infeasible
t_lower = 0;
t_upper = 1;
while sol.problem == 1
t_upper = t_upper*2;
sol = solvesdp([C>=0, A <= t_upper*B, B>=0],[],ops);
end
end
%else
% error('Hmm, you''d better robustify the code!')
%end

% OK, on a les limites , on commence la bissection
while t_upper-t_lower >= 1e-3
t_mid = (t_upper+t_lower)/2;
sol = solvesdp([C>=0, A <= t_mid*B, B>=0],[],ops);
if sol.problem == 0
t_upper = t_mid;
else
t_lower = t_mid;
end
%disp([t_lower t_upper]);
end
*************************************************************
I have obtained the following results :
iter seconds digits c*x b*y
21 0.4 Inf -7.8725074213e+001 -7.8725071780e+001
|Ax-b| = 4.7e-010, [Ay-c]_+ = 8.5E-010, |x|= 1.9e+004, |y|= 7.0e+001

Detailed timing (sec)
Pre IPM Post
3.120E-002 3.900E-001 0.000E+000
Max-norms: ||b||=1, ||c|| = 7.989655e+001,
Cholesky |add|=0, |skip| = 0, ||L.L|| = 2.44718.
ans =
yalmiptime: 0.225999999999999
solvertime: 0.326999999999998
info: 'Successfully solved (SeDuMi-1.1)'
problem: 0
dimacs: [NaN NaN NaN NaN NaN NaN]
X =
8.938487227006524 3.738119536210800
3.738119536210800 46.878546847633686
Pi =
9.480842159422313 5.051249067025251
5.051249067025251 69.244229620159260
ans =
22.444410959621575
Fi =
-0.573452926049265
ans =
-1
1
ans =
-0.286726463024632 + 2.742185592942308i
-0.286726463024632 - 2.742185592942308i

but the exact results (in the article SOF stabilization : an ILMI approach - Cao et al. , 1998) are :
ai =
-0.0377
Pi =
9.2557 6.48
6.4888 92.9376
Fi =
-0.7369

How is wrong with my code , please ?