Question: fsolve cannot solve this equations and return them.

hello

this is my program and fsolve for low intensity solve the equations but for high intensity cannot solve why?

this is my code:

ep0 := 1/(4*3.14);

el := 8.54*10^(-2);

hbar := 1;

vf := 1/300;

kb := 1;

tem := 2.586*10^(-2);

ci := 1;

p := 1.458*10^16;

beta := 2;

ai := 7.1*10^(-4);

bi := ai/sqrt(3);

enph := .196;

d := enph/(kb*tem);

n0 := 1/(exp(enph/(kb*tem))-1);

gama := hbar*vf;

intensity:=9000000

w := 7.28;

impurity := 7.2*10^3;

g := hbar*beta/(bi^2*sqrt(2*p*enph));

aa := g^2*(n0+1)/(2*Pi*hbar*gama^2);

bb := g^2*n0/(2*Pi*hbar*gama^2);

cc := 2/(Pi*gama^2);

l := (1*hbar)*w/(2*kb*tem);

u := el^2*intensity/(32*w*hbar^2);

[fsolve({op([((enph*ln(1+exp(c+enph/(kb*tem)))/(kb*tem)-polylog(2, -exp(c))+polylog(2, -exp(c+enph/(kb*tem))))*enph*(kb*tem)^2-(enph^2*ln(1+exp(c+enph/(kb*tem)))/(kb^2*tem^2)+2*enph*polylog(2, -exp(c+enph/(kb*tem)))/(kb*tem)+2*polylog(3, -exp(c))-2*polylog(3, -exp(c+enph/(kb*tem))))*(kb*tem)^3+(-exp(b)*enph*ln(1+exp(c+enph/(kb*tem)))+exp(c+d)*enph*ln(1+exp(b-d+enph/(kb*tem)))+exp(b)*kb*tem*polylog(2, -exp(c))-exp(c+d)*kb*tem*polylog(2, -exp(b-d))-exp(b)*kb*tem*polylog(2, -exp(c+enph/(kb*tem)))+exp(c+d)*kb*tem*polylog(2, -exp(b-d+enph/(kb*tem))))*enph*(kb*tem)^2/((exp(b)-exp(c+d))*kb*tem)+(exp(b)*enph^2*ln(1+exp(c+enph/(kb*tem)))-exp(c+d)*enph^2*ln(1+exp(b-d+enph/(kb*tem)))+2*exp(b)*enph*kb*tem*polylog(2, -exp(c+enph/(kb*tem)))-2*exp(c+d)*enph*kb*tem*polylog(2, -exp(b-d+enph/(kb*tem)))+2*exp(b)*kb^2*tem^2*polylog(3, -exp(c))-2*exp(c+d)*kb^2*tem^2*polylog(3, -exp(b-d))-2*exp(b)*kb^2*tem^2*polylog(3, -exp(c+enph/(kb*tem)))+2*exp(c+d)*kb^2*tem^2*polylog(3, -exp(b-d+enph/(kb*tem))))*(kb*tem)^3/((exp(b)-exp(c+d))*kb^2*tem^2))*bb+u*(1/(1+exp(-l-c))-1/((1+exp(-l-c))*(1+exp(l-b))))-(((1*enph)*(enph-2*kb*tem*ln(1+exp(-b+enph/(kb*tem))))/(2*kb^2*tem^2)+2*kb^2*tem^2*(-polylog(2, -exp(-b+enph/(kb*tem)))+polylog(2, -cosh(b)+sinh(b))))*enph*(kb*tem)^2-(enph^2*(enph-3*kb*tem*ln(1+exp(-b+enph/(kb*tem))))-6*kb^2*tem^2*(enph*polylog(2, -exp(-b+enph/(kb*tem)))+kb*tem*(-polylog(3, -exp(-b+enph/(kb*tem)))+polylog(3, -cosh(b)+sinh(b)))))*(kb*tem)^3/(3*kb^3*tem^3)-(-exp(b)*enph^2+exp(c+d)*enph^2-2*exp(c+d)*enph*kb*tem*ln(1+exp(-b+enph/(kb*tem)))+2*exp(b)*enph*kb*tem*ln(1+exp(-c-d+enph/(kb*tem)))+2*exp(c+d)*kb^2*tem^2*polylog(2, -exp(-b))-2*exp(b)*kb^2*tem^2*polylog(2, -exp(-c-d))-2*exp(c+d)*kb^2*tem^2*polylog(2, -exp(-b+enph/(kb*tem)))+2*exp(b)*kb^2*tem^2*polylog(2, -exp(-c-d+enph/(kb*tem))))*enph*(kb*tem)^2/((2*(-exp(b)+exp(c+d)))*kb^2*tem^2)-(exp(b)*enph^3-exp(c+d)*enph^3+3*exp(c+d)*enph^2*kb*tem*ln(1+exp(-b+enph/(kb*tem)))-3*exp(b)*enph^2*kb*tem*ln(1+exp(-c-d+enph/(kb*tem)))+6*exp(c+d)*enph*kb^2*tem^2*polylog(2, -exp(-b+enph/(kb*tem)))-6*exp(b)*enph*kb^2*tem^2*polylog(2, -exp(-c-d+enph/(kb*tem)))+6*exp(c+d)*kb^3*tem^3*polylog(3, -exp(-b))-6*exp(b)*kb^3*tem^3*polylog(3, -exp(-c-d))-6*exp(c+d)*kb^3*tem^3*polylog(3, -exp(-b+enph/(kb*tem)))+6*exp(b)*kb^3*tem^3*polylog(3, -exp(-c-d+enph/(kb*tem))))*(kb*tem)^3/((3*(-exp(b)+exp(c+d)))*kb^3*tem^3))*aa-u*(1/(1+exp(l-b))-1/((1+exp(-l-c))*(1+exp(l-b)))) = 0, -cc*polylog(2, -exp(b))+cc*polylog(2, -exp(-c))-impurity = 0])}, {op([b, c])})];

thank you.