MATLAB: Solving system of integral equations – Error FSOLVE cannot continue

fsolvesystem of integral equations

Dear all,
I am trying to solve a system of complicated integral equations in MATLAB.
The equations and initial values are attached in the figure.
Indeed, I want to solve Eqs. (35) and (37) with considering Eq. (36) vanishes in vacuum. Of course the next step will be solving all three equations with together.
I am using Fsolve routine in Matlab R2015a to obtaine "m" and "Delta_MF", But some thing is wrong and I get error massage in the command window. I would really appreciate that if someone could help me to solve my equations.
T=50.;
mu=0.;
mustar = 0.;
Lambda = 602.3;
% Lambda is the upper limit of integral.
GsLambda2 = 2.319;
Gs = 0.00000639;
m0 = 5.5;
Gv = 0.5 * Gs;
Gd = 0.75 * Gs;
mgap = @(b) integral(@(p) 4.0/pi.^2*Gs*b(1)*1.0/(sqrt(p.^2+b(1).^2))*((1.-2.*((exp ((sqrt ((sqrt(p.^2+b(1).^2) - mustar).^2 + (abs (b(2))).^2))/T) + 1).^-1))*(sqrt(p.^2+b(1).^2) - mustar)/(sqrt ((sqrt(p.^2+b(1).^2) - mustar).^2 + (abs (b(2))).^2))+(1-2.*((exp ((sqrt ((sqrt(p.^2+b(1).^2) + mustar).^2 + (abs (b(2))).^2))/T) + 1).^-1))*(sqrt(p.^2+b(1).^2) + mustar)/(sqrt ((sqrt(p.^2+b(1).^2) + mustar).^2 + (abs (b(2))).^2))+((exp (-(sqrt(p.^2+b(1).^2) + mustar)/T) + 1).^-1)-((exp ((sqrt(p.^2+b(1).^2) - mustar)/T) + 1).^-1)), 0, 602.3) - (b(1)-m0);
Deltagap = @(b) integral(@(p) 4.0/pi.^2*Gd*b(2)*((1-2.*((exp ((sqrt ((sqrt(p.^2+b(1).^2) - mustar).^2 + (abs (b(2))).^2))/T) + 1).^-1))/(sqrt ((sqrt(p.^2+b(1).^2) - mustar).^2 + (abs (b(2))).^2))+(1-2.*((exp ((sqrt ((sqrt(p.^2+b(1).^2) + mustar).^2 + (abs (b(2))).^2))/T) + 1).^-1))/(sqrt ((sqrt(p.^2+b(1).^2) + mustar).^2 + (abs (b(2))).^2))), 0, 602.3) - b(2);
B = fsolve(@(b) [mgap(b), Deltagap(b)], rand(2,1))
Error using  / 
Matrix dimensions must agree.
Error in
@(p)4.0/pi.^2*Gs*b(1)*1.0/(sqrt(p.^2+b(1).^2))*((1.-2.*((exp((sqrt((sqrt(p.^2+b(1).^2)-mustar).^2+(abs(b(2))).^2))/T)+1).^-1))*(sqrt(p.^2+b(1).^2)-mustar)/(sqrt((sqrt(p.^2+b(1).^2)-mustar).^2+(abs(b(2))).^2))+(1-2.*((exp((sqrt((sqrt(p.^2+b(1).^2)+mustar).^2+(abs(b(2))).^2))/T)+1).^-1))*(sqrt(p.^2+b(1).^2)+mustar)/(sqrt((sqrt(p.^2+b(1).^2)+mustar).^2+(abs(b(2))).^2))+((exp(-(sqrt(p.^2+b(1).^2)+mustar)/T)+1).^-1)-((exp((sqrt(p.^2+b(1).^2)-mustar)/T)+1).^-1))
Error in integralCalc/iterateScalarValued (line 314)
                fx = FUN(t);
Error in integralCalc/vadapt (line 132)
            [q,errbnd] = iterateScalarValued(u,tinterval,pathlen);
Error in integralCalc (line 75)
        [q,errbnd] = vadapt(@AtoBInvTransform,interval);
Error in integral (line 88)
Q = integralCalc(fun,a,b,opstruct);
Error in
@(b)integral(@(p)4.0/pi.^2*Gs*b(1)*1.0/(sqrt(p.^2+b(1).^2))*((1.-2.*((exp((sqrt((sqrt(p.^2+b(1).^2)-mustar).^2+(abs(b(2))).^2))/T)+1).^-1))*(sqrt(p.^2+b(1).^2)-mustar)/(sqrt((sqrt(p.^2+b(1).^2)-mustar).^2+(abs(b(2))).^2))+(1-2.*((exp((sqrt((sqrt(p.^2+b(1).^2)+mustar).^2+(abs(b(2))).^2))/T)+1).^-1))*(sqrt(p.^2+b(1).^2)+mustar)/(sqrt((sqrt(p.^2+b(1).^2)+mustar).^2+(abs(b(2))).^2))+((exp(-(sqrt(p.^2+b(1).^2)+mustar)/T)+1).^-1)-((exp((sqrt(p.^2+b(1).^2)-mustar)/T)+1).^-1)),0,602.3)-(b(1)-m0)
Error in @(b)[mgap(b),Deltagap(b)]
Error in fsolve (line 219)
            fuser = feval(funfcn{3},x,varargin{:});
Error in gapequations (line 41)
B = fsolve(@(b) [mgap(b), Deltagap(b)], rand(2,1))
Caused by:
    Failure in initial user-supplied objective function evaluation. FSOLVE cannot continue.

Best Answer

Failure in initial user-supplied objective function evaluation. FSOLVE cannot continue.
Before you use FSOLVE, you must make sure that your objective function code works. 
T=50.;
mu=0.;
mustar = 0.;
Lambda = 602.3;
% Lambda is the upper limit of integral.
GsLambda2 = 2.319;
Gs = 0.00000639;
m0 = 5.5;
Gv = 0.5 * Gs;
Gd = 0.75 * Gs;
mgap = @(b) integral(@(p) 4.0/pi.^2*Gs*b(1)*1.0/(sqrt(p.^2+b(1).^2))*((1.-2.*((exp ((sqrt ((sqrt(p.^2+b(1).^2) - mustar).^2 + (abs (b(2))).^2))/T) + 1).^-1))*(sqrt(p.^2+b(1).^2) - mustar)/(sqrt ((sqrt(p.^2+b(1).^2) - mustar).^2 + (abs (b(2))).^2))+(1-2.*((exp ((sqrt ((sqrt(p.^2+b(1).^2) + mustar).^2 + (abs (b(2))).^2))/T) + 1).^-1))*(sqrt(p.^2+b(1).^2) + mustar)/(sqrt ((sqrt(p.^2+b(1).^2) + mustar).^2 + (abs (b(2))).^2))+((exp (-(sqrt(p.^2+b(1).^2) + mustar)/T) + 1).^-1)-((exp ((sqrt(p.^2+b(1).^2) - mustar)/T) + 1).^-1)), 0, 602.3) - (b(1)-m0);
Deltagap = @(b) integral(@(p) 4.0/pi.^2*Gd*b(2)*((1-2.*((exp ((sqrt ((sqrt(p.^2+b(1).^2) - mustar).^2 + (abs (b(2))).^2))/T) + 1).^-1))/(sqrt ((sqrt(p.^2+b(1).^2) - mustar).^2 + (abs (b(2))).^2))+(1-2.*((exp ((sqrt ((sqrt(p.^2+b(1).^2) + mustar).^2 + (abs (b(2))).^2))/T) + 1).^-1))/(sqrt ((sqrt(p.^2+b(1).^2) + mustar).^2 + (abs (b(2))).^2))), 0, 602.3) - b(2);
fun=@(b) [mgap(b), Deltagap(b)];
fun(rand(2,1))
Error using /
Matrix dimensions must agree.

Error in solution (line 11)
mgap = @(b) integral(@(p) 4.0/pi.^2*Gs*b(1)*1.0/(sqrt(p.^2+b(1).^2))*((1.-2.*((exp ((sqrt ((sqrt(p.^2+b(1).^2) - mustar).^2 + (abs (b(2))).^2))/T) + 1).^-1))*(sqrt(p.^2+b(1).^2) - mustar)/(sqrt ((sqrt(p.^2+b(1).^2) - mustar).^2 + (abs (b(2))).^2))+(1-2.*((exp ((sqrt ((sqrt(p.^2+b(1).^2) + mustar).^2 + (abs (b(2))).^2))/T) + 1).^-1))*(sqrt(p.^2+b(1).^2) + mustar)/(sqrt ((sqrt(p.^2+b(1).^2) + mustar).^2 + (abs (b(2))).^2))+((exp (-(sqrt(p.^2+b(1).^2) + mustar)/T) + 1).^-1)-((exp ((sqrt(p.^2+b(1).^2) - mustar)/T) + 1).^-1)), 0, 602.3) - (b(1)-m0);

Error in integralCalc/iterateScalarValued (line 314)
                fx = FUN(t);

Error in integralCalc/vadapt (line 132)
            [q,errbnd] = iterateScalarValued(u,tinterval,pathlen);

Error in integralCalc (line 75)
        [q,errbnd] = vadapt(@AtoBInvTransform,interval);

Error in integral (line 88)
Q = integralCalc(fun,a,b,opstruct);

Error in solution (line 11)
mgap = @(b) integral(@(p) 4.0/pi.^2*Gs*b(1)*1.0/(sqrt(p.^2+b(1).^2))*((1.-2.*((exp ((sqrt ((sqrt(p.^2+b(1).^2) - mustar).^2 + (abs (b(2))).^2))/T) + 1).^-1))*(sqrt(p.^2+b(1).^2) - mustar)/(sqrt ((sqrt(p.^2+b(1).^2) - mustar).^2 + (abs (b(2))).^2))+(1-2.*((exp ((sqrt ((sqrt(p.^2+b(1).^2) + mustar).^2 + (abs (b(2))).^2))/T) + 1).^-1))*(sqrt(p.^2+b(1).^2) + mustar)/(sqrt ((sqrt(p.^2+b(1).^2) + mustar).^2 + (abs (b(2))).^2))+((exp (-(sqrt(p.^2+b(1).^2) + mustar)/T) + 1).^-1)-((exp ((sqrt(p.^2+b(1).^2) - mustar)/T) + 1).^-1)), 0, 602.3) - (b(1)-m0);

Error in solution (line 14)
fun=@(b) [mgap(b), Deltagap(b)];
Related Question