please help me, i wrote the following code. but all results is NaN and give me empty diagram . where the problem in my code or where is the error ?

% [fxxdc]=therioticalfunction(1,0.2)

function [ fxxdc]=therioticalfunction(ld,a)

[w1,w2]=meshgrid(0:pi/20:pi,0:pi/20:pi);

fxxdc=abs((((ld*exp(ld)*(3*ld + exp(2*ld) - 2*exp(ld) - 3*ld*exp(ld) + ld^2 + ld^2*exp(ld) + 1))/(exp(ld) - 1)^3 + (ld*exp(ld)*((exp(-ld) - 1)*a - 2*ld*exp(-ld)*a^2)*(ld - exp(ld) + 1))/((exp(ld) - 1)^2*((1 - exp(-ld))*a^2 + (exp(-ld) - 1)*a)))*((a^3*exp(i*w1)*exp(-i*w2)*(exp(-ld) - 1)^2)/((a*exp(-i*w2)*(exp(-ld) - 1) + 1)*(exp(i*w1)*(exp(-ld) - 1)*a^2 + 1)) + (a^3*exp(i*w2)*exp(-i*w1)*(exp(-ld) - 1)^2)/((a*exp(-i*w1)*(exp(-ld) - 1) + 1)*(exp(i*w2)*(exp(-ld) - 1)*a^2 + 1)) + (a^3*exp(i*(w1 + w2))*exp(i*w1)*(exp(-ld) - 1)^2)/((a*exp(i*(w1 + w2))*(exp(-ld) - 1) + 1)*(exp(i*w1)*(exp(-ld) - 1)*a^2 + 1)) + (a^3*exp(i*(w1 + w2))*exp(i*w2)*(exp(-ld) - 1)^2)/((a*exp(i*(w1 + w2))*(exp(-ld) - 1) + 1)*(exp(i*w2)*(exp(-ld) - 1)*a^2 + 1)) + (a^3*exp(-i*(w1 + w2))*exp(-i*w1)*(exp(-ld) - 1)^2)/((a*exp(-i*w1)*(exp(-ld) - 1) + 1)*(exp(-i*(w1 + w2))*(exp(-ld) - 1)*a^2 + 1)) + (a^3*exp(-i*(w1 + w2))*exp(-i*w2)*(exp(-ld) - 1)^2)/((a*exp(-i*w2)*(exp(-ld) - 1) + 1)*(exp(-i*(w1 + w2))*(exp(-ld) - 1)*a^2 + 1))) - ((ld*exp(ld)*(3*ld + exp(2*ld) - 2*exp(ld) - 3*ld*exp(ld) + ld^2 + ld^2*exp(ld) + 1))/(exp(ld) - 1)^3 + (ld*exp(ld)*((exp(-ld) - 1)*a - 2*ld*exp(-ld)*a^2)*(ld - exp(ld) + 1))/((exp(ld) - 1)^2*((1 - exp(-ld))*a^2 + (exp(-ld) - 1)*a)))*((a^2*exp(-i*(w1 + w2))*(exp(-ld) - 1))/(exp(-i*(w1 + w2))*(exp(-ld) - 1)*a^2 + 1) + (a^2*exp(i*w1)*(exp(-ld) - 1))/(exp(i*w1)*(exp(-ld) - 1)*a^2 + 1) + (a^2*exp(i*w2)*(exp(-ld) - 1))/(exp(i*w2)*(exp(-ld) - 1)*a^2 + 1)) + (ld*exp(ld)*(3*ld + exp(2*ld) - 2*exp(ld) - 3*ld*exp(ld) + ld^2 + ld^2*exp(ld) + 1))/(exp(ld) - 1)^3 - (ld*exp(ld)*((a*exp(i*(w1 + w2))*(exp(-ld) - 1))/(a*exp(i*(w1 + w2))*(exp(-ld) - 1) + 1) + (a*exp(-i*w1)*(exp(-ld) - 1))/(a*exp(-i*w1)*(exp(-ld) - 1) + 1) + (a*exp(-i*w2)*(exp(-ld) - 1))/(a*exp(-i*w2)*(exp(-ld) - 1) + 1))*(3*ld + exp(2*ld) - 2*exp(ld) - 3*ld*exp(ld) + ld^2 + ld^2*exp(ld) + 1))/(exp(ld) - 1)^3 - (ld*exp(ld)*((exp(-ld) - 1)*a - 2*ld*exp(-ld)*a^2)*(ld - exp(ld) + 1)*((a^2*exp(i*(w1 + w2))*exp(i*w1)*(exp(-ld) - 1)^2)/((a*exp(i*(w1 + w2))*(exp(-ld) - 1) + 1)*(a*exp(i*w1)*(exp(-ld) - 1) + 1)) + (a^2*exp(i*(w1 + w2))*exp(i*w2)*(exp(-ld) - 1)^2)/((a*exp(i*(w1 + w2))*(exp(-ld) - 1) + 1)*(a*exp(i*w2)*(exp(-ld) - 1) + 1)) + (a^2*exp(-i*(w1 + w2))*exp(-i*w1)*(exp(-ld) - 1)^2)/((a*exp(-i*(w1 + w2))*(exp(-ld) - 1) + 1)*(a*exp(-i*w1)*(exp(-ld) - 1) + 1)) + (a^2*exp(-i*(w1 + w2))*exp(-i*w2)*(exp(-ld) - 1)^2)/((a*exp(-i*(w1 + w2))*(exp(-ld) - 1) + 1)*(a*exp(-i*w2)*(exp(-ld) - 1) + 1)) + (a^2*exp(i*w1)*exp(-i*w2)*(exp(-ld) - 1)^2)/((a*exp(i*w1)*(exp(-ld) - 1) + 1)*(a*exp(-i*w2)*(exp(-ld) - 1) + 1)) + (a^2*exp(i*w2)*exp(-i*w1)*(exp(-ld) - 1)^2)/((a*exp(i*w2)*(exp(-ld) - 1) + 1)*(a*exp(-i*w1)*(exp(-ld) - 1) + 1))))/((exp(ld) - 1)^2*((1 - exp(-ld))*a^2 + (exp(-ld) - 1)*a)) + (ld*exp(ld)*((exp(-ld) - 1)*a - 2*ld*exp(-ld)*a^2)*((a*exp(-i*(w1 + w2))*(exp(-ld) - 1))/(a*exp(-i*(w1 + w2))*(exp(-ld) - 1) + 1) + (a*exp(i*w1)*(exp(-ld) - 1))/(a*exp(i*w1)*(exp(-ld) - 1) + 1) + (a*exp(i*w2)*(exp(-ld) - 1))/(a*exp(i*w2)*(exp(-ld) - 1) + 1))*(ld - exp(ld) + 1))/((exp(ld) - 1)^2*((1 - exp(-ld))*a^2 + (exp(-ld) - 1)*a)))/(4*pi^2));

mesh(w1,w2,fxxdc)

end

[Rafael Hernandez-Walls]

I think your code can improve. You can try changing all multiplications (*) to (. *) And all divisions (/) to (./) and you will see that your code improves. This because you are working with two matrices (w1 and w2).

[mohammed elmenshawy]

Thank you very much , Really correct answer.