@mmcdara 

Many thanks

@adel-00 

lambda:=0.1:N:=5:M:=sqrt(N(N+1)):omegap:=10:phi:=0:

var:={n(t),u(t)}:
dsys:={diff(n(t),t)=-2*(n(t)-N)+(u(t)-abs(M)*exp(I*phi))*exp(-2*I*omegap*t/lambda)+conjugate((u(t)-abs(M)*exp(I*phi))*exp(-2*I*omegap*t/lambda)),diff(u(t),t)=-2*(1-I*delta)*u(t)+2*(n(t)-N)*exp(-2*I*omegap*t/lambda)+2*abs(M)*exp(I*phi)}:
res1:=dsolve(dsys union {n(0)=0,u(0)=0},numeric,output=listprocedure,parameters=[delta]):

    
WV := eval(n(t), res1):

tlo,thi := 0.0, 2.0;
U0lo,U0hi := -10, 10;(m,n) := 51,51;
    
str := time[real]():
M:=Matrix(m,n,datatype=float[8]):
for j from 1 to n do
  u0 := U0lo+(j-1)*(U0hi-U0lo)/(n-1);
  WV(parameters=[delta=u0]);
  for i from 1 to m do
    T := tlo+(i-1)*(thi-tlo)/(m-1);
    try
      M[i,j] := WV(T);
    catch:
      # mostly maxfun complaint for t above some value.
      M[i,j] := Float(undefined);
    end try;
  end do:
end do:
plots:-surfdata(M, tlo..thi, U0lo..U0hi,
                labels=["t","delta","n(t)"]);

Found similar to my problem  which explained by Carl

@mmcdara 

thanks

want to plot 3d as (n(t,delta))

@Carl Love 

Thanks very much...

@Carl Love 

very usefull

u := proc (x, t) options operator, arrow; gamma1*tanh(2*Pi-t-x) end proc;
v := proc (x, t) options operator, arrow; gamma2*tanh(2*Pi+t+x) end proc;
A := 1/4*(-2*u(x, t)*(diff(v(x, t), x))+2*beta1*u(x, t)*(diff(u(x, t), t, t))+beta2*u(x, t)^2*(u(x, t)^2+v(x, t)^2)-2*alpha1*u(x, t)*(diff(v(x, t), t, t, t))-alpha2*u(x, t)*(u(x, t)^2+v(x, t)^2)*(diff(v(x, t), t))+2*v(x, t)*(diff(u(x, t), x))+2*beta1*v(x, t)*(diff(v(x, t), t, t))+beta2*v(x, t)^2*(u(x, t)^2+v(x, t)^2)+2*alpha1*v(x, t)*(diff(u(x, t), t, t, t))+alpha2*v(x, t)*(u(x, t)^2+v(x, t)^2)*(diff(u(x, t), t))):
            u := (x, t) -> gamma1 tanh(2 Pi - t - x)
            v := (x, t) -> gamma2 tanh(2 Pi + t + x)
Coeffs:= proc(e::algebraic, V::{list,set}(name))
local T, C;
    try C:= coeffs(e, V, T)
    catch "invalid argument":
        C:= coeffs(collect(e, V, ':-distributed'), V, T)
    end try;
    table([T]=~[C])
end proc
:
V:= [gamma1,gamma2,alpha1, alpha2, beta1, beta2]:
C:= Coeffs(A, V):
r1:= add(v*int(C[v], [x= -3..3, t= -10..10], numeric), v= [indices](C, nolist)):

eq1:=(diff(r1,gamma1)):eq2:=(diff(r1,gamma2)):
sol = evalf(solve({eq1, eq2}, {gamma1, gamma2}));

Root_Z.. appear in the sol

Thanks for ue comments

still i got error

Download Carl_code.mwCarl_code.mw

@Carl Love 

@Carl Love 

one more thing 

if 

u := proc (x, t) options operator, arrow; a*tanh(2*Pi-t-x) end proc;
v := proc (x, t) options operator, arrow; b* tanh(2*Pi+t+x) end proc;

two more constants are added

@Carl Love 

Thanks Carl

is 

r1:= add(v*int(C[v], [x= -10..10, t= -3..3], numeric), v= V);

double integral??

u := proc (x, t) options operator, arrow; tanh(2*Pi-t-x) end proc;
v := proc (x, t) options operator, arrow; tanh(2*Pi+t+x) end proc;
A := 1/4*(-2*u(x, t)*(diff(v(x, t), x))+2*beta1*u(x, t)*(diff(u(x, t), t, t))+beta2*u(x, t)^2*(u(x, t)^2+v(x, t)^2)-2*alpha1*u(x, t)*(diff(v(x, t), t, t, t))-alpha2*u(x, t)*(u(x, t)^2+v(x, t)^2)*(diff(v(x, t), t))+2*v(x, t)*(diff(u(x, t), x))+2*beta1*v(x, t)*(diff(v(x, t), t, t))+beta2*v(x, t)^2*(u(x, t)^2+v(x, t)^2)+2*alpha1*v(x, t)*(diff(u(x, t), t, t, t))+alpha2*v(x, t)*(u(x, t)^2+v(x, t)^2)*(diff(u(x, t), t)));
f := (x,t) -> A:
r1:=int(int(A,x=-10..10),t=-3..3,numeric):

it takes time 

@acer 

u := proc (x, t) options operator, arrow; tanh(2*Pi-t-x) end proc;
v := proc (x, t) options operator, arrow; tanh(2*Pi+t+x) end proc;
A := a*u(x, t)*(diff(u(x, t), t))-u(x, t)*(diff(u(x, t), x, x, x))+b*v(x, t)*(diff(v(x, t), t))-v(x, t)*(diff(v(x, t), x, x, x));
f := (x,t) -> A:
r1:=evalf(int(int(A,x=-10..10),t=-10..10,numeric)):

@Polovodov 

thanks..

if ( a) is inside the integrals and cant go outside!! is there a way to eval the integral without assigning (a)..

for other example i have

@acer 

Thanks for your responds

first: y is a complex; lets say y=a+ib, So y^2 is Re(y)^2+Im(y)^2

the isuue is all integrations contains deltao(unknown), so how can I plot y^2 against deltao 

is there  a way to do it without loop?

it is not at each point of the vectror yr is assigned to sls

for example for 0.3 in vertical axis should be verticle coulmn (200 )

@Carl Love 

Thanks

But why in matlab produce the figure different?

what shall i add in maple in the  plot option to get the same result..

@Ramakrishnan 

Thanks

xa and xb are functions of t

P1:=logplot([xa,xb],0..1,axes=boxed,title=tit,color=black,font=[2,3,18],thickness=2,tickmarks=[3,2],titlefont=[SYMBOL,14],font=[1,1,18],linestyle=[1,3,2]);

but still the same issue not like the matlab plot