Question: plot a function divided by its maximum

Hi,

1-Triying to plot a function divided by its maximum value,sometimes it works with some parameters that means, the max.value of the plot is 1.

But when i change the data the max. value in the plot in graeter than 1 which is wrong!! should be 1.

dont know why??

2- Changing different data in the parameters, the programme takes long long time then i stop it?

please help me with these two problems.

restart:
>
------------------------- Defining the nature of the variables used ----------------------
assume(T,real):Digits:=25:n:=1:tau:=Pi:
theta:=0:phi:=0:
lambda:=n;Omega:=1:Gamma:=0.01:
--------------------- Input---------------------------------
1

J1

term1:=(exp((Gamma+I*d)*tau)-1)/(2*(Gamma+I*d)):
Ak1:=d->(exp((Gamma+I*(d+0.5*Omega-2*lambda*k/Gamma))*tau)-1)/(Gamma+I*(d+0.5*Omega-2*lambda*k/Gamma))+(exp((Gamma+I*(d-0.5*Omega+2*lambda*k/Gamma))*tau)-1)/(Gamma+I*(d-0.5*Omega+2*lambda*k/Gamma)):
Ak2:=d->(exp((Gamma+I*(d-0.5*Omega-2*lambda*k/Gamma))*tau)-1)/(Gamma+I*(d-0.5*Omega-2*lambda*k/Gamma))+(exp((Gamma+I*(d+0.5*Omega+2*lambda*k/Gamma))*tau)-1)/(Gamma+I*(d+0.5*Omega+2*lambda*k/Gamma)):
term2:=(evalf(-0.25*sum(BesselJ(k,Omega*Gamma/(4*n))*Ak1(d)+BesselJ(k,-(Omega*Gamma)/(4*n))*Ak2(d),k=0..50))):
J1:=(term1+term2):
J1mod:=(Re(J1))^2+(Im(J1))^2:
###### J2#########################
Ak1:=d->(exp((Gamma+I*(d+0.5*Omega-2*lambda*k/Gamma))*tau)-1)/(Gamma+I*(d+0.5*Omega-2*lambda*k/Gamma))-(exp((Gamma+I*(d-0.5*Omega+2*lambda*k/Gamma))*tau)-1)/(Gamma+I*(d-0.5*Omega+2*lambda*k/Gamma)):
Ak2:=d->(exp((Gamma+I*(d-0.5*Omega-2*lambda*k/Gamma))*tau)-1)/(Gamma+I*(d-0.5*Omega-2*lambda*k/Gamma))-(exp((Gamma+I*(d+0.5*Omega+2*lambda*k/Gamma))*tau)-1)/(Gamma+I*(d+0.5*Omega+2*lambda*k/Gamma)):

J2:=(evalf(-0.25*sum(BesselJ(k,Omega*Gamma/(4*n))*Ak1(d)+BesselJ(k,-Omega*Gamma/(4*n))*Ak2(d),k=0..100))):
######################

J2mod:=(Re(J2))^2+(Im(J2))^2:
J3 same as J1differ in sign
term1:=(exp((Gamma+I*d)*tau)-1)/(2*(Gamma+I*d)):
Ak1:=d->(exp((Gamma+I*(d+0.5*Omega-2*lambda*k/Gamma))*tau)-1)/(Gamma+I*(d+0.5*Omega-2*lambda*k/Gamma))+(exp((Gamma+I*(d-0.5*Omega+2*lambda*k/Gamma))*tau)-1)/(Gamma+I*(d-0.5*Omega+2*lambda*k/Gamma)):
Ak2:=d->(exp((Gamma+I*(d-0.5*Omega-2*lambda*k/Gamma))*tau)-1)/(Gamma+I*(d-0.5*Omega-2*lambda*k/Gamma))+(exp((Gamma+I*(d+0.5*Omega+2*lambda*k/Gamma))*tau)-1)/(Gamma+I*(d+0.5*Omega+2*lambda*k/Gamma)):
term2:=(evalf(0.25*sum(BesselJ(k,Omega*Gamma/(4*n))*Ak1(d)+BesselJ(k,-Omega*Gamma/(4*n))*Ak2(d),k=0..100))):
J3:=term1+term2:
J3mod:=(Re(J3))^2+(Im(J3))^2:
J4 same as J2 but -0.25-->2

J4:=-2*J2:
######################

J4mod:=(Re(J4))^2+(Im(J4))^2:

calculate the spectrum

Spec:=d->(exp(-2*Gamma*tau)*(J1mod*cos(theta/2)^2+J2mod+J3mod*sin(theta/2)^2-0.5*Re(J3*J4*sin(theta)*exp(I*phi))+0.5*Re(J1*J4*sin(theta)*exp(-I*phi)))):

with(plots):

tit:=sprintf("l=%g,W=%g,G=%g",lambda,Omega,Gamma):
Smax1:=max(seq(evalf(Spec(d)),d=-100..100)):
plot(evalf(Spec(d)/Smax1),d=-15..15,axes=boxed,title=tit,color=black,font=[2,3,18],thickness=2,tickmarks=[3,3],titlefont=[SYMBOL,14],font=[1,1,18],linestyle=1);