@acer 

Hi 

I found error in the integration so i did correct the expressions but I got computation interrupted:

here the code below plz I appreciated if you can give me comments

restart: # d is real (plot of Spec agianst d)
    
kernelopts(version):    

term1:=(exp((1+I*d)*Gamma*tau)*(1-cos(2*Omega1))+cos(2*Omega1)*exp((1+I*d)*Gamma*t0)-exp(-(1+I*d)*Gamma*t0))/(2*(1+I*d)):
term2:=-1/2*evalf(Int(exp((1+I*d)*Gamma*x)*cos(Omega1*(1+erf(x))),x=-3..3)):
J1:=evalf(term1+term2):

J1mod:=(Re(J1))^2+(Im(J1))^2:

###### J2#########################
A2:=Int(exp((1+I*d)*Gamma*x)*sin(Omega1*(1+erf(x))),x=-3..3):
A3:=sin(2*Omega1)*(exp((1+I*d)*Gamma*tau)-exp((1+I*d)*Gamma*t0))/(1+I*d):
J2:=-I*(A2+A3):
J2mod:=(Re(J2))^2+(Im(J2))^2:
#J3 same as J1differ in sign
term1:=(exp((1+I*d)*Gamma*tau)*(1+cos(2*Omega1))-cos(2*Omega1)*exp((1+I*d)*Gamma*t0)+exp(-(1+I*d)*Gamma*t0))/(2*(1+I*d)):
term2:=0.5*int(exp((1+I*d)*Gamma*x)*cos(Omega1*(1+erf(x))),x=-3..3):
J3:=term1+term2:
J3mod:=(Re(J3))^2+(Im(J3))^2:
J4:=-J2:
J4mod:=(Re(J4))^2+(Im(J4))^2:

#Spec:=J1mod+J2mod:
Spec:=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):
tau:=Pi:tau0:=1:Omega1:=10;Gamma:=1:t0:=3:theta:=Pi/3:phi:=Pi/3:
                          Omega1 := 10
tit:=sprintf("W=%g,G=%g,t=%g,q=%g,f=%g",Omega1,Gamma,tau,theta,phi);
         tit := "W=10,G=1,t=3.14159,q=1.0472,f=1.0472"
P1:=plot(Spec,d=-350..350,axes=boxed,numpoints=200,title=tit,color=black,
         font=[2,3,18],thickness=2,tickmarks=[3,3],gridlines=false,
         labels=["",""],
         titlefont=[SYMBOL,14],font=[1,1,18],linestyle=1);
Warning,  computation interrupted
Normalize:= proc(P::specfunc(anything, PLOT))
  local A,Smax1;
  A:= op([1,1], P);
  Smax1:= max(A[..,2]);
  if A::list then A:= Matrix(A) end if;
  A[..,2]:= A[..,2]/Smax1;
  subsop([1,1]= A, P);
end proc:
P1:= Normalize(P1):
for kk from 2 to 3 do
  tau:= 0.2*kk*Pi;
  P||kk:= plot(Spec,d=-350..350,axes=boxed,numpoints=200,title=tit,color=black,
               font=[2,3,18],thickness=2,tickmarks=[3,3],gridlines=false,
               labels=["",""],
               titlefont=[SYMBOL,14],font=[1,1,18],linestyle=1);
  P||kk:= plottools:-translate(Normalize(P||kk), 0, kk-1)
od:
display([P||(1..3)],view=[-10..10,0..3]);

@acer 

it is really run quicker many thanks

@acer 

              "Normalised spectrum for all cases"
restart:
------------------------- Defining the nature of the variables used ----------------------
assume(theta,real):assume(phi,real):assume(d,real):assume(tau,real):

                           tau0 := 1
J1
term1:=(exp((1+I*d)*Gamma*tau)*(1-cos(2*Omega1))+cos(2*Omega)*exp((1+I*d)*Gamma*t0)-exp(-(1+I*d)*Gamma*t0))/(2*(1+I*d)):

term2:=-0.5*evalf(int(exp((1+I*d)*Gamma*x)*cos(Omega1*(1+erf(tau))),x=-1..1)):
J1:=(term1+term2):
J1mod:=(Re(J1))^2+(Im(J1))^2:

###### J2#########################
A2:=evalf(int(exp((1+I*d)*Gamma*x)*sin(Omega1*(1+erf(tau))),x=-1..1)):
A3:=sin(2*Omega1)*(exp((1+I*d)*Gamma*tau)-exp((1+I*d)*Gamma*t0))/(1+I*d):

J2:=-I*(A2+A3):
######################

J2mod:=(Re(J2))^2+(Im(J2))^2:

J3 same as J1differ in sign
term1:=(exp((1+I*d)*Gamma*tau)*(1+cos(2*Omega1))-cos(2*Omega)*exp((1+I*d)*Gamma*t0)+exp(-(1+I*d)*Gamma*t0))/(2*(1+I*d)):

term2:=0.5*evalf(int(exp((1+I*d)*Gamma*x)*cos(Omega1*(1+erf(tau))),x=-1..1)):

J3:=term1+term2:
J3mod:=(Re(J3))^2+(Im(J3))^2:
J4 =- J2 
J4:=-J2:
J4mod:=(Re(J4))^2+(Im(J4))^2:
calculate the spectrum
Spec:= unapply((J1mod+J2mod,d)):
#Spec:= unapply((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))),d):
with(plots):

Omega:=1:tau:=Pi:tau0:=1:Omega1:=0.5*Omega*tau0*sqrt(Pi):Gamma:=0.05:t0:=0.75:theta:=0:phi:=0:
#tit:=sprintf("l=%g,W=%g,G=%g,q=%g,f=%g",lambda,Omega,Gamma,theta,Phi):
P1:=plot((Spec),-350..350,axes=boxed,numpoints=200,title=tit,color=black,font=[2,3,18],thickness=2,tickmarks=[3,3],titlefont=[SYMBOL,14],font=[1,1,18],linestyle=1);

Normalize:= proc(P::specfunc(anything, PLOT))
local A,Smax1;
A:= op([1,1], P);
Smax1:= max(A[..,2]);
if A::list then A:= Matrix(A) end if;
A[..,2]:= A[..,2]/Smax1;
subsop([1,1]= A, P)
end proc:

P1:= Normalize(P1):
for kk from 2 to 5 do
tau0:= kk*Pi;
P||kk:= plot(Spec,-350..350,axes=boxed,numpoints=200,title=tit,color=black,font=[2,3,18],thickness=2,tickmarks=[3,3],titlefont=[SYMBOL,14],font=[1,1,18],linestyle=1);
P||kk:= plottools:-translate(Normalize(P||kk), 0, kk-1)
od:

display([P||(1..5)],view=[-350..350,0..5]);

that is the code that the integration in the program.. but it isn;t work

@acer 

hi if i make changes where

evalf(int(exp((5+I*d)*x)*sin(1+erf(x)),x=-1.2..1.2));

as d is areal number

in my opinoin I try to:

sin(erf) is equal to 0.5( exp(erf)-iexp(-erf))where erf(x)= 2/sqrt(pi) sum(x^2n+1)n!(2n+1)

so if I find the term exp(ax+bx^3) will be brilliant 

@Christian Wolinski 

to find integration of sin(erf(t)) e^(a+ib)t I need explicit formula

many thanks of ur response

@Carl Love I am not sure if I can exapnd it either with Bessel function or somthing simmilar

if i integrate sin(erf(t)) e^(a+ib)t

@acer thanks I got it

@Preben Alsholm 

and with that

restart:Digits:=70:
------------------------- Defining the nature of the variables used ----------------------
N:=0:M:=0:N1:=1+N:w:=10:

ini1:= x(0)=0.5,y(0)=0.5,z(0)=0;
            ini1 := x(0) = 0.5, y(0) = 0.5, z(0) = 0
var:={x(t),y(t),z(t)}: 
dsys:={diff(z(t),t)=-(N1+M*cos(2*w*t))*z(t)-1+f*(x(t)+y(t)), diff(x(t),t)=-(N1-I*w-2*M*exp(-2*I*w*t))*x(t)-f*(N1+(z(t)))-2*f*M*exp(2*I*w*t),diff(y(t),t)=-(N1+I*w-2*M*exp(2*I*w*t))*y(t)-f*(N1+(z(t)))-2*f*M*exp(-2*I*w*t)}:
zd:=subs(dsys,diff(z(t),t));
               zd := -z(t) - 1 + f (x(t) + y(t))
res:=dsolve(dsys union {x(0)=0.5,y(0)=0.5,z(0)=0},numeric,output=listprocedure):
Warning, The use of global variables in numerical ODE problems is deprecated, and will be removed in a future release. Use the 'parameters' argument instead (see ?dsolve,numeric,parameters)
#tit:=sprintf("F=%g,N=%g",f,N):
plot3d(z(t,f), f=-1..1, t=0..3, grid=[29,49]);

@acer 

Many thanks for ur response 

If I want to plot z(t) 

CodeTools:-Usage( plot3d(zfun(t,f), f=-1..1, t=0..3, grid=[29,49]) );

@Preben Alsholm 

N:=0:M:=0:N1:=1+N:w:=10:

ini1:= x(0)=0.5,y(0)=0.5,z(0)=0;
            ini1 := x(0) = 0.5, y(0) = 0.5, z(0) = 0
var:={x(t),y(t),z(t)}: 
dsys:={diff(z(t),t)=-(N1+M*cos(2*w*t))*z(t)-1+f*(x(t)+y(t)), diff(x(t),t)=-(N1-I*w-2*M*exp(-2*I*w*t))*x(t)-f*(N1+(z(t)))-2*f*M*exp(2*I*w*t),diff(y(t),t)=-(N1+I*w-2*M*exp(2*I*w*t))*y(t)-f*(N1+(z(t)))-2*f*M*exp(-2*I*w*t)}:
zd:=subs(dsys,diff(z(t),t));
                 zd := -z(t) - 1 + x(t) + y(t)
res:=dsolve(dsys union {x(0)=0.5,y(0)=0.5,z(0)=0},numeric,output=listprocedure):

I mean plot 3d of z(t),t,f 

many thanks

and for plotting 3d g(z(t),t,f)

@Carl Love 

Dear Carl

If want to plot

f:=epsilon=(lambda1+sqrt(-lambda2^2+Y^2*(lambda3^2+(lambda4-lambda5*Y^2)^2)));

epsilon against the sqrt(Y)

is it correct like this

r1:=plot([rhs(f), sqrt(Y), sqrt(Y)= 0..1],axes=boxed,thickness=2,color=black,font=[1,1,20],tickmarks=[3, 3],linestyle=1,labels= [typeset(epsiloni), sqrt(Y)]);

@Carl Love 

restart:with(plots):
Digits:=35:
------------------------- Defining the nature of the variables used ----------------------
assume(Y,real);
Omega:=5*Pi:
gamma1:=8*Pi:
gamma2:=0.002*Pi:
x1:=100*Pi:
omega2:=200*Pi:
gamma0:=0.2*Pi:
#Delta:=25*Pi:
G:=20*Pi:

#tit:=sprintf("G=%g,delta1=%g,gamma1=%g",G,delta1,gamma1):
lambda1:=(1/(2*Pi))^2*(G*Omega*gamma0/(2*(0.25*gamma0^2+Delta^2))):
lambda2:=(1/(2*Pi))^2*(-G*Omega*Delta/((0.25*gamma0^2+Delta^2))):
lambda3:=(1/(2*Pi))^2*gamma1+lambda1:
lambda4:=(1/(2*Pi))^2*(0.2*omega2-G^2*Delta/(0.25*gamma0^2+Delta^2)):
lambda5:=(1/(2*Pi))^2*(2*x1^2*omega2/((omega2^2+gamma2^2))):
epsilon-(lambda1+sqrt(-lambda2^2+Y*(lambda3^2+(lambda4-lambda5*Y)^2)))^2:
#implicitplot(f,epsilon=0..200,Y=0..1,numpoints=1000,axes=boxed,thickness=2,color=black,font=[1,1,20],tickmarks=[3, 3],linestyle=1):
implicitplot3d(f, epsilon=0..20,Y = 0 .. 0.4,Delta=-25*Pi..25*Pi,labels=[E1,Y,Delta],tickmarks=[3,3,3]);
with(Student[Calculus1]):contourplot(f,Delta =-25*Pi .. 25*Pi,Y=0..1,contours=[0],axes=boxed,thickness=2,color=black,font=[1,1,18],tickmarks=[5, 5],linestyle=1,view=[0..1,0..0.6]);

@acer 

sorry

I thought it is not clear here

i have done the implicit plot for 3d, but i dont know how to do the contour for the implicit expression