@Axel Vogt 

R:=evalc(1/(x^3+x^2+x+2));plot(R)

complex plot what i mean

@mmcdara 

Thanks for your efforts.. I tried to follow your comments but it takes time to run  the loop. Thanks again

@Carl Love 

Hi Carl,

Sorry for not clicking on the thumbs up. your  valuble comments are helpful

How can I make normaliztion ( as all values are nearly zero for large d)

I.e. how to divide the function (that we  plot) on its max

@Carl Love 

No reason for computing over -20.. 20 it just a try.. the benifit is to run faster.

Thanks,

I find a way also in loop as:

restart:with(plots):
Omega:=0.01:no:=1:Delta:=0:tau0:=1:c:=0:gamma1:=0:j:=1:k1:=0:t:=10:xr:=1:Gamma:=0.1:

for i from -80 to 80 do
d[i]:=i:
b:=-gamma1*tau0+I*tau0*Delta-2*a*(k-1)*xr:
a:=1+I*c;
c1:=sqrt(conjugate(a))/tau0:
c2:=0.5*((conjugate(b)/sqrt(conjugate(a)))):
lambda1[i]:=2*(Gamma-I*d[i])*(t-k1)/j+(Gamma-I*d[i])^2:
lambda2[i]:=t*(1-j)*k1/j+1/sqrt(2)*(Gamma-I*d[i]):
lambda3[i]:=c1*(t-k1)/j+c1*(Gamma-I*d[i])+c2:
J1:=sqrt(Pi)/sqrt(2)*(1-erf(lambda2[i])):
J1mod:=(Re(J1))^2+(Im(J1))^2:
g1:=0.5*sqrt(Pi)*tau0*exp(c2^2)*exp(-conjugate(a)*((k-1)*xr)^2)/(sqrt(conjugate(a))):
F2[i]:=evalf(-sqrt(2)*Int(exp(-x^2)*erf(sqrt(2)*c1*x+lambda3[i]),x=-lambda2[i]..5)):
J2:=sum(g1*J1+g1*F2[i],k=1..1):
J2mod:=(((Re(J2))^2+(Im(J2))^2)):
f[d[i]]:=J2mod:

end:

ptsN1 := [seq([d[i], f[d[i]]], i = -20 .. 20)]:
P11 := plot(ptsN1, color = black, linestyle = solid, thickness = 2, color = black, axes = boxed, labelfont = ["HELVETICA", 14], labels = ["d", "S_M(d)"])

which gives the same plot as yours...

I guess the problem with the singularity as u menshioned.

Thanks for ur valuable comments

1- the lower integral is complex.

2- sum can’t change it to eval because it is from k=1 to 1 (for know ) later will increase  the value of k.

3- I’ll  change infinity to large value as 100.

Now it takes time to plot.. any suggestions plz