Why do I keep getting "no solution found" message?

Hello. I am trying to solve a system of two equations that one of them is an algebraic equation and the second one is integral equation, but changing the initial guesses I keep getting "no solution found" or when the equation is solved, the answer is imaginary that can not be. for both unknown variables, only positive values are acceptable. I would really appreciate that if someone can help me in this.

here is the solver with all input parameters in the equations:

mpi = 138.0;

fpi = 93.0;

msigma = 370.58;

m0 = 790.0;

g1 = 13.0;

g2 = 6.96;

gp = g1+g2;

gm = g1-g2;

gomega = 6.80;

momega = 783;

gro = 7.98;

epsilon = (mpi.^2*fpi);

L2 = 0.5*(msigma.^2-3.*mpi.^2.);

L4 = (msigma.^2-mpi^2)/(2.*fpi.^2.);

f=zeros(150,7);

for j = 1:1:150

nb = j.*0.01;

omega = (gomega/(momega.^2)).*nb.*(197.3).^3;

%

sigmagap = @(S) [2./(3.*pi.^2).*((S(1).^2 - 1./4.*(sqrt(gp.^2.*S(2).^2 + 4.*m0.^2) + gm.*S(2)).^2) + (S(1).^2 - 1./4.*(sqrt(gp.^2.*S(2).^2 + 4.*m0.^2) - gm.*S(2)).^2))-(nb.*197.3^3);

integral(@(p) 1.0/(2.*pi.^2)*p.^2.* (((gp + ((S(2).*gm.^2)./sqrt ((gp.^2.*S(2).^2)+4.*790.^2))).* ...

((gm.*S(2) + sqrt ((gp.^2.*S(2).^2)+4.*790.^2))./( sqrt (p.^2.+ (1./4.*(sqrt(gp.^2*S(2).^2+4.*790.^2)+gm.*S(2)).^2.).^2 ))))), 0, sqrt(S(1).^2 - 1./4.*(sqrt(gp.^2.*S(2).^2 + 4.*m0.^2) + gm.*S(2)).^2), 'ArrayValued',1) + ...

integral(@(p) 1.0/(2.*pi.^2)*p.^2.* (( (-gm + ((S(2).*gp.^2)./sqrt ((gp.^2.*S(2).^2)+4.*790.^2))).* ...

((-gm.*S(2) + sqrt ((gp.^2.*S(2).^2)+4.*790.^2))./ ...

( sqrt (p.^2.+ (1./4.*(sqrt(gp.^2*S(2).^2+4.*790.^2)-gm.*S(2)).^2.).^2 ))))), 0, sqrt(S(1).^2 - 1./4.*(sqrt(gp.^2.*S(2).^2 + 4.*m0.^2) - gm.*S(2)).^2), 'ArrayValued',1) - epsilon - L2.*S(2) + L4.*S(2).^3];

S0 = [900; 90];

options = optimoptions('fsolve','Display','iter');

[S,fval] = fsolve(sigmagap,S0,options);

mstarp = 0.5*(sqrt(gp.^2.*S(2).^2.+4.*m0.^2)+gm.*S(2));

mstarm = 0.5*(sqrt(gp.^2.*S(2).^2.+4.*m0.^2)-gm.*S(2));

mub = S(1)+gomega.*omega;

f(j,1) = nb;

f(j,2) = S(2);

f(j,3) = omega;

f(j,4) = mstarp;

f(j,5) = mstarm;

f(j,6) = S(1);

f(j,7) = mub;

end

Answer 1

[Torsten]

Maybe a better code structure to debug the values "fsolve" gets from "sigmagap".

You think that you know that "sigmagap" returns complex values, don't you ?

mpi = 138.0;
fpi = 93.0;
msigma = 370.58;
m0 = 790.0;
g1 = 13.0;
g2 = 6.96;
gp = g1+g2;
gm = g1-g2;
gomega = 6.80;
momega = 783; 
gro = 7.98;
epsilon = (mpi.^2*fpi);
L2 = 0.5*(msigma.^2-3.*mpi.^2.);
L4 = (msigma.^2-mpi^2)/(2.*fpi.^2.);
f=zeros(150,7);
S0 = [900;90];
for j = 1:1:1
    nb = j.*0.01;
    sigmagap(S0,gp,m0,gm,nb,epsilon,L2,L4)
    %options = optimset('Display','none'); 
    [S,fval] = fsolve(@(S)sigmagap(S,gp,m0,gm,nb,epsilon,L2,L4),S0);%,options)
    mstarp = 0.5*(sqrt(gp.^2.*S(2).^2.+4.*m0.^2)+gm.*S(2));
    mstarm = 0.5*(sqrt(gp.^2.*S(2).^2.+4.*m0.^2)-gm.*S(2));
    omega = (gomega/(momega.^2)).*nb.*(197.3).^3;
    mub = S(1)+gomega.*omega;
    f(j,1) = nb;
    f(j,2) = S(2);
    f(j,3) = omega;
    f(j,4) = mstarp;
    f(j,5) = mstarm;
    f(j,6) = S(1);
    f(j,7) = mub;
    S0 = S;
end

ans = 
   1.0e+05 *

  -1.7066 + 0.0000i  -3.9500 - 7.6862i

No solution found.

fsolve stopped because the relative size of the current step is less than the
value of the step size tolerance squared, but the vector of function values
is not near zero as measured by the value of the function tolerance.

function res = sigmagap(S,gp,m0,gm,nb,epsilon,L2,L4)
  res(1) = 2./(3.*pi.^2).*((S(1).^2 - 1./4.*(sqrt(gp.^2.*S(2).^2 + 4.*m0.^2) + gm.*S(2)).^2) + ...
           (S(1).^2 - 1./4.*(sqrt(gp.^2.*S(2).^2 + 4.*m0.^2) - gm.*S(2)).^2))-(nb.*197.3^3);
  res(2) =  integral(@(p) 1.0/(2.*pi.^2)*p.^2.* ...
            (((gp + ((S(2).*gm.^2)./sqrt ((gp.^2.*S(2).^2)+4.*790.^2))).* ...
            ((gm.*S(2) + sqrt ((gp.^2.*S(2).^2)+4.*790.^2))./ ...
            ( sqrt (p.^2.+ (1./4.*(sqrt(gp.^2*S(2).^2+4.*790.^2)+gm.*S(2)).^2.).^2 ))))), 0, ...
            sqrt(S(1).^2 - 1./4.*(sqrt(gp.^2.*S(2).^2 + 4.*m0.^2) + gm.*S(2)).^2), 'ArrayValued',1) + ...
            integral(@(p) 1.0/(2.*pi.^2)*p.^2.* ...
            (( (-gm + ((S(2).*gp.^2)./sqrt ((gp.^2.*S(2).^2)+4.*790.^2))).* ...
            ((-gm.*S(2) + sqrt ((gp.^2.*S(2).^2)+4.*790.^2))./ ...
            ( sqrt (p.^2.+ (1./4.*(sqrt(gp.^2*S(2).^2+4.*790.^2)-gm.*S(2)).^2.).^2 ))))), 0, ...
            sqrt(S(1).^2 - 1./4.*(sqrt(gp.^2.*S(2).^2 + 4.*m0.^2) - gm.*S(2)).^2), 'ArrayValued',1) ...
            - epsilon - L2.*S(2) + L4.*S(2).^3;
end