The following error occurred converting from function_handle to double:
Error using double
Conversion to double from function_handle is not possible.

clear all
nn=3;
mm=3;
totalatom=mm*nn;
vf=1e6;
el=1.6e-19;
energy=1e-19;
hbar=1.054e-34;
omega=energy/hbar;
ep0=8.854e-12;
refindex=2.71+complex(1.41);
permitrel=5.356+complex(0,1)*7.642;
permit=ep0*permitrel;
c=3e8;
intensity=10e4;
power=1;
pulsedur=1e-15;
photonflux=intensity/energy;
photonnumber=power*pulsedur/energy;
a0=(hbar*refindex*photonflux/(2*photonnumber*omega*permit*c));
% Hamiltonian Of Graphene Without Perturbation
h=ones(totalatom,totalatom);
for i=1:totalatom;
    for j=1:totalatom;
        if j==i-1;
            h(i,j)=2;
        elseif j==i+1;
            h(i,j)=2;
       % elseif j==i+3;
        %    h(i,j)=2;
        else
            h(i,j)=0;
        end
    end
end


%syms e;
dive=1e-8;
green=@(e) e*double(eye(totalatom))-h-dive*double(eye(totalatom))*complex(0,1);
gret=@(e) eye(totalatom)/green;
gless=@(e)complex(0,1)*gret; % Remark: Here The Fermi-Dirac Distribution Is Not Write It Will Product To Finall Answer
gless=@(e) subs(gless,e,e-energy);

sigma=ones(totalatom,totalatom);
coupling=0;
for l=1:totalatom
    for m=1:totalatom
        for p=1:totalatom
            if p==l+1
                mlp=complex(0,1)*el*vf*a0*square(3)/2;
            elseif p==l-1
                mlp=-complex(0,1)*el*vf*a0*square(3)/2;
            elseif p==l
                continue
            elseif p==l+2
                break
            end
                for q=1:totalatom
                    if q==m+1
                        mqm=complex(0,1)*el*vf*a0*square(3)/2;
                    elseif q==m-1
                        mqm=-complex(0,1)*el*vf*a0*square(3)/2;
                    elseif q==m
                        mqm=0;
                    elseif q==m+2
                        break
                    end
                    coupling=@(e) mlp*mqm*gless(p,q)+coupling;
                end
        end
        sigma(l,m)=@(e) coupling*photonnumber;
        coupling=0;
    end
end
%Lesser Green Function Photon
glessphoton=@(e) gret*sigma*ctranspose(gret);
diag=@(e) trace(glessphoton);
resp=integral(@(e) diag,0,inf);