Self-consistent solution of integral equations using fsolve

% numerical self-consistent calculation 
clc
clear variables
warning('off')
fileID = fopen('data.txt','w');
KbT = logspace(-4,2,21);    
for i = 1:length(KbT)
Dd = 10.0;
tn = 0.2;
ed = -0.5;
delta = 10^(-6);
a = KbT(i)./tn;
syms g Gr11 Ga11 G11r G11a
x = (1:10000000);
x(1)=0.4; %initial guess
n = 1;
while true     
    m =  100; % number of data points for integration wrt to 'z' using trapezoidal rule
    z=linspace(-10.0,10.0,m);
    y = zeros(1,100);
    
    for k=1:m      
        Gr = fun(z(k),Dd,ed,tn,KbT(i),delta,x(n));
        %y = (-1./pi).*((1./2).*(1-tanh(z(k)./(2.*KbT(i))))).*imag(Gr11);
        y = (-1./pi).*(1./(exp(z(k)./KbT(i))+1)).*imag(Gr);
        Wanted_sol(k) = double(y);     
    end
    %x(n+1) = integral(@(t) interp1(z,Wanted_sol,t,'linear','extrap'), z(1), z(end),'ArrayValued',true);
     x(n+1) = quadgk(@(t) interp1(z,Wanted_sol,t,'linear','extrap'), z(1), z(end),'RelTol',0,'AbsTol',1e-9);
    if (abs(x(n+1)-x(n))<0.001)
       ndown = x(n+1);
       nup = ndown;
       m =  10; % number of data points for integration wrt to 'z' using simpson rule
       omega=linspace(-5.*KbT(i),5.*KbT(i),m);
       f1 = zeros(1,10);
       f2 = zeros(1,10);
       f3 = zeros(1,10);
        for j = 1:length(omega)
           Gr = fun(omega(j),Dd,ed,tn,KbT(i),delta,nup);
           Ga = conj(Gr);
           T1 =  (tn.^2).*(Gr.*Ga);
           fdd = (1./KbT(i)).*(exp(omega(j)./KbT(i))./(exp(omega(j)./KbT(i))+1).^2);
           f1(j) =  fdd.*T1;
        end
       % Use quadgk to integrate the data
        L01 =  2.*quadgk(@(t) interp1(omega,double(f1),t,'linear','extrap'), omega(1), omega(end),'RelTol',0,'AbsTol',1e-9);
       % Use Gaussian quadrature to integrate the data
       %L01 =  2.*integral(@(t) interp1(omega,f1,t,'linear','extrap'), omega(1), omega(end),'ArrayValued',true);
       G =  L01;
      fprintf(fileID,'%8.6e %8.6e %8.6e\n',[a,G,nup]');
      fprintf('%8.6e %8.6e %8.6e\n',[a,G,nup]') 
       break
    end
    x(n)=x(n+1);
    n=n+1;
end
end
fclose(fileID);
%setting the Matlab figure
d = load('data.txt');
a = d(:,1);
b = d(:,2);
c = d(:,3);
semilogx(a,b,'o-K','Linewidth', 2.0,'Markersize',4.0)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function Gr = fun(z,Dd,ed,tn,KbT,delta,x)
options = optimoptions('fsolve','Display','off','TolFun',1e-9,'TolX',1e-9);
% Integration limit
    lower_lmt = -10.0;
    upper_lmt = 10.0;
    y_0 = [0.1; 0.1];
    % self-consistent equations
    F = @(y) double([
    y(1)-(tn./pi).*quadgk(@(z1) ((((1./(exp(z1./KbT)+1))).*(((1-x-(y(1)))./(z1-ed+(tn./pi).*log(abs((Dd-z1)./(Dd+z1)))-1i.*tn+1i.*tn.*(y(1))-(y(2))))))...
    ./(z-z1+1i.*delta.*heaviside(z)-1i.*delta.*heaviside(-z))), lower_lmt, upper_lmt, 'RelTol',0,'AbsTol',1e-9);
    y(2)-(tn./pi).*quadgk(@(z1) (((1./(exp(z1./KbT)+1)).*(1+1i.*tn.*(((1-x-(y(1)))./(z1-ed+(tn./pi).*log(abs((Dd-z1)./(Dd+z1)))-1i.*tn+1i.*tn.*(y(1))-(y(2)))))))...
    ./(z-z1+1i.*delta.*heaviside(z)-1i.*delta.*heaviside(-z))), lower_lmt, upper_lmt, 'RelTol',0,'AbsTol',1e-9)...
    ]);
    sol = fsolve(F,y_0,options);
    eng_1 = vpa(sol(1));
    eng_2 = vpa(sol(2));
    g0 =   z-ed+(tn./pi).*log(abs((Dd-z)./(Dd+z)))+1i.*tn;
    Gr =  ((1-x-eng_1)./(g0-eng_2-1i.*tn.*eng_1));
end