You brought up a similar question here:

http://www.mapleprimes.com/questions/207342-The-Sample-Option-In-Odeplot-

@GeorgeReynolds It seems that you misunderstood my latest reply.
My example was intended to show that an error (a strange one) results when using D-format, i.e. when pdeD is used, whereas pde (in diff-format) gives an answer.

@leosh The variable _Z appears in the RootOf expressions in answer

indets([answer],RootOf);
     {RootOf(-L^2*cos(xi)^2+_Z^2+w^2), RootOf(_Z^2-2), RootOf(-L^2+_Z^2+w^2), RootOf(-L^2+_Z^2+2*w^2)}

RootOf(_Z^2-2) (as an example) stands for the roots of the polynomial _Z^2-2, thus stands for sqrt(2) and -sqrt(2).

Use allvalues to get all the roots:

map(allvalues,{answer});

You can avoid the RooOfs by using
answer := solve({equ1, equ2, equ3, equ4, equ5, equ6}, {alpha, sigma, xi, lambda__1, lambda__2, beta__f},explicit=true);

Notice that the two variable arctan is used:  arctan(y,x)
In the help page for arctan you find this:
"For real arguments x, y, the two-argument function arctan(y, x), computes the principal value of the argument of the complex number x + I y , so -Pi < arctan(y, x) and arctan(y, x) <= Pi. "

I had problems with equ6. When I copied your input line to a 1D math line I got the following weird stuff (well, I hate 2D math input with a passion):
`#mi("equ6")` := `#mi("L")`*`#mi("sin",fontstyle = "normal")`(sigma)*cos(xi)-w*alpha*v^2*sin(sigma)/(g*l*cos(xi)^2)

After having corrected that I could concentrate on the problem.
The answer is that solve simply doesn't find any solution, therefore it returns NULL (i.e. nothing).

Could you edit this so we only get the input? As it is now it is unreadable.

In Maple go to Edit/Remove Output/From Worksheet

After that copy all the lines and paste them into the MaplePrimes editor.

@Rouben Rostamian  And the logarithmic plot:

plots:-logplot(abs(myx1(t) - myx2(t)), t=0..4*Pi);

@GeorgeReynolds 
1. Consider the following simple example:
restart;
pde:=diff(u(x,t),t)=diff(u(x,t),x,x);
pdsolve(pde,{u(x,0)=x,u(0,t)=0,D[1](u)(1,t)=1},numeric); ## OK
pdeD:=convert(pde,D);
pdsolve(pdeD,{u(x,0)=x,u(0,t)=0,D[1](u)(1,t)=1},numeric);

Error, (in pdsolve/numeric/process_PDEs) number of dependent variables and number of PDE must be the same

This is clearly a limitation in the implementation of pdsolve/numeric. It is stated explicitly in the help:
?pdsolve,numeric
"The PDE system PDEsys must be specified in diff notation."

There is no complaint from the symbolic solver:
pdsolve(pdeD); #Separates variables
pdsolve({pdeD,u(x,0)=x,u(0,t)=0,D[1](u)(1,t)=1}); #No error, but no result (as I expected)

2. No, surely not. I'm pretty certain that it uses discretization in space as well as in time.

@GeorgeReynolds You cannot use numerical solution unless the parameters are given numerical values. You have 4 parameters with no values given. One is I which you surely don't mean to be the imaginary unit that Maple takes it to be.
There are several solutions to this problem. The simplest (to me) is to rename your I to II or whatever you like.
There are 3 more:

indets(PDE1,name) minus {xi,vartheta};

Output:      {E, a, omega}

########### After that you will run into other problems. To avoid one of them correct your line
convert(PDE1, diff);
## to
PDE1 := convert(PDE1, diff);
## Next problem: Maple's numerical solver is time based. Let us assume that vartheta is your time. Thus it starts from initial values given on the xi interval xi=0..1 (in your case) and proceeds in time from there satisfying whatever boundary conditions are given at xi=0 and xi=1. You don't have any initial values, only boundary conditions.

##Example, where I have replaced I with II:
PDE1 := convert(PDE1, diff);
params:=indets(PDE1,name) minus {xi,vartheta};
parvals:=params=~1; #Setting all 4 parameters to 1
eval(PDE1, parvals); #Just to have a look
#ic:={V(xi,0)=xi,D[2](V)(xi,0)=0}; #An example
ic:={V(xi,0)=xi^2*(1-xi),D[2](V)(xi,0)=0}; #Maybe a better example
pds := pdsolve(eval(PDE1, parvals), bc union ic,numeric);
pds:-plot(vartheta=1);

@GeorgeReynolds To use Worksheet interface and 1D math input (in the session or globally) go to (in Windows):


Tools/Options/Display/Input display. Choose Maple Notation.  (i.e. 1D math).
Tools/Options/Interface/Default format for new workswheets. Choose Worksheet.

Finally press either Apply to Session or Apply Globally.

The changes only affects new worksheets.
There is no danger in choosing the latter since the change can always be reverted.  

What I see is (for the initial value problem, as an example):

y(x) = JacobiSN((1/2)*sqrt(2)*x+InverseJacobiSN(1, I), I)

and after a % I get

y(x) = JacobiSN((1/2)*sqrt(2)*x+EllipticK(I), I)

How did you get those sn's in the first place? Is there in fact something missing in your worksheet?

@Markiyan Hirnyk Since you need a restart in between, it is difficult to compare two exact result that are so big:

But just try

restart;
int(......);
evalf(%);

restart;
int(......);
evalf(%);

There you surely don't have numerical integration, so the cause is roundoff.

@Konstantin@ Yes, if int cannot do the integration (as is the case here), then it uses numerical integration.

Try

int(exp(sin(x)+x+exp(x)),x=0..1);

@Bendesarts The code I used:

restart;
K:=Matrix([<0, -1, 1, -1>,<-1, 0, -1, 1>,<-1, 1, 0,-1>,<1, -1, -1,0>]);
omega[sw]:=beta/(1-beta)*omega[s];
for i to 4
do
r[i]:=sqrt((u[i](t)/(L/2))^2+(v[i](t)/H)^2):
omega[i]:=omega[st]/(1+exp(b*v[i](t)))+omega[sw]/(1+exp(-b*v[i](t))):
Equ[i]:=diff(u[i](t),t)=Au*(1-r[i]^2)*u[i](t)+omega[i]*(L/2)/H*v[i](t):
Eqv[i]:=diff(v[i](t),t)=Av*(1-r[i]^2)*v[i](t)+omega[i]*(L/2)/H*v[i](t)+(K.<seq(v[i](t),i=1..4)>)[i]:
EqSys[i]:=[Equ[i],Eqv[i]]:
end do;
paramsGeo:=L=0.015,H=0.015,beta=0.5,Vf=0.3;
omegaS:=eval(Pi*Vf/L, [paramsGeo]);
paramsCycle:=omega[s]=omegaS,Au=1,Av=1,b=100;
params:=paramsGeo,paramsCycle;

sys:=map(op,eval({seq(EqSys[i],i=1..4)},{params}));
#ic:={u[1](0)=0.8, v[1](0)=0,u[2](0)=0.8, v[2](0)=0,u[3](0)=0.8, v[3](0)=0,u[4](0)=0.8, v[4](0)=0};
ic:={u[1](0)=0.8, v[1](0)=0.1,u[2](0)=0.8, v[2](0)=0.1,u[3](0)=0.8, v[3](0)=0.1,u[4](0)=0.8, v[4](0)=0.1};
res:=dsolve(eval(sys,omega[st]=50) union ic,numeric);
plots:-odeplot(res,[u[1](t),v[1](t)],0..1);
plots:-odeplot(res,[seq([t,v[i](t)+i/10],i=1..4)],0..1,thickness=2);
#The last plot:

I don't see any text in the body of the question. How come Carl can see it?

@GambiaMan What is it you want plotted: Only the points corresponding to the times in your sample? If that is the case your option numpoints=2000 confuses me.