Get rid of eps which will do nothing useful.
Then
expression y;
y=0;
for j=1:2
for i=1:2
y=y+a(i,j).*log( (s(i,j)*fai) /(a(i,j)+eps) );
end
end
maximize(y)
can be replaced with
maximize(sum(sum(-rel_entr(a,a+s*fai))))
The constraints can be more simply written as
sum(s,1) <= min(1,L/(B * e_ele);
a >= s
0 <= a <= 1
0 <= s <= 1
Follow the instructions at CVXQUAD: How to use CVXQUAD's Pade Approximant instead of CVX's unreliable Successive Approximation for GP mode, log, exp, entr, rel_entr, kl_div, log_det, det_rootn, exponential cone. CVXQUAD's Quantum (Matrix) Entropy & Matrix Log related functions . Use Mosek 9.x with CVX 2.2 if available to you, otherwise, install CVXQUAD with its exponential.m replacement. No changes to your code are necessary because CVXQUAD processes rel_entr as is.