To see the limit cycle mentioned by me under the discussion using DEplot you don't need all these initial points.

Using one suffices. I chose one not too far away from the 3rd equilibrium point:

sol:=solve({eq1=0,eq2=0,eq3=0,eq4=0,eq5=0},{x,y,z,v,w});
pts:=map(subs,[sol],[x,y,z,v,w]);
eps:=.1:
ini:=[x,y,z,v,w](0)=~pts[3]+[eps$5];
res:=dsolve({Eq1,Eq2,Eq3,Eq4,Eq5} union {ini[]},numeric,maxfun=0);
plots:-odeplot(res,[x(t),y(t)],500..511);
XYZ:=plots:-odeplot(res,[x(t),y(t),z(t)],500..511):
YZV:=plots:-odeplot(res,[y(t),z(t),v(t)],500..511):
ZVW:=plots:-odeplot(res,[z(t),v(t),w(t)],500..511):
plots:-display(Array([XYZ,YZV,ZVW]));

@J4James I don't know what your objective is, but often people are interested in the long term behavior of the solutions (example: the Lorenz strange attractor).

If that is also your interest then you should use a time range like t=t1..t2, where t1 is somewhat large and t2-t1 is large enough to give an idea of the long term behavior. In case of a limit cycle that time difference should be just large enough to represent the cycle. In case of a more complicated long term behavior it is not so clear.

But start by trying  t = 250..260 or similar. While experimenting you might as well use stepsize=0.1.
Replace the default thickness with the one suggested by Carl.

There is graphical evidence of a limit cycle. To add more excitement I included z in the plot and used t=500..510.
Notice the need for increasing maxfun (actually I set it at infinity by maxfun=0):

DEtools[DEplot3d]([Eq1,Eq2,Eq3,Eq4,Eq5],[x(t),y(t),z(t),v(t),w(t)],t=500..510,initialset,stepsize=0.1,color=blue,linecolor=magenta,axes=boxed,scene=[x,y,z],x=5.2..7.2,y=0.2..1.6,z=0..0.15,obsrange=false,thickness=0,maxfun=0);

@J4James I now remember that there is an option ,obsrange=false.

Use that and keep your original range.

@Carl Love I think you are right. I tried splitting the initial set in two according to the following procedure Q.
Notice that the range for x has been widened to begin with 5.

Q:=proc(ini) local x1,y1,a,b,c,d;
    a:=5; b:=7.2; c:=0.2; d:=1.6;
    x1,y1:=op(subs(ini,[x(0),y(0)])); (x1-a)*(x1-b)<=0 and (y1-c)*(y1-d)<=0
end proc;
Q(initialset[166]); #Test

iniset1,iniset2:=selectremove(Q,initialset):
nops~([iniset1,iniset2]);

DEtools[DEplot]([Eq1,Eq2,Eq3,Eq4,Eq5],[x(t),y(t),z(t),v(t),w(t)],t=0..30,x=5..7.2,y=0.2..1.6,iniset1 ,stepsize=0.1,color=blue,linecolor=magenta,axes=boxed,scene=[x,y]);

DEtools[DEplot]([Eq1,Eq2,Eq3,Eq4,Eq5],[x(t),y(t),z(t),v(t),w(t)],t=0..30,x=5..7.2,y=0.2..1.6,iniset2 ,stepsize=0.1,color=blue,linecolor=magenta,axes=boxed,scene=[x,y]);

@J4James I tried the following, where I only took the last 11 of the initial conditions and only used time 0..30 (my patience is limited too). I removed restrictions on x and y and also the arrows option as it doesn't make much sense here (and is ignored anyway).

It worked fine.

DEtools[DEplot]([Eq1,Eq2,Eq3,Eq4,Eq5],[x(t),y(t),z(t),v(t),w(t)],t=0..30,initialset[170..180],stepsize=0.01,color=blue,linecolor=magenta,axes=boxed,scene=[x,y]);

@Vasudeva I see your point.

So if we let alpha = diff(beta(t),t) we get

ode:=diff(alpha(beta),beta)=1/alpha(beta)*sin(beta)*(r*alpha(beta)^2+g)/(2*r*cos(beta)-J);

#Using dsolve with the 'implicit' option, we get

dsolve(ode,implicit);

#and without:
sol:=dsolve(ode);

#Thus we have for the first one

ode2:=subs({alpha(beta)=diff(beta(t),t),beta=beta(t)},sol[1]);
dsolve(ode2);
res:=value(%);
odetest(res,eq);
#The second:
ode2b:=subs({alpha(beta)=diff(beta(t),t),beta=beta(t)},sol[2]);
dsolve(ode2b);
resb:=value(%);
odetest(resb,eq);

Without trying to understand what is going on I notice that you are using equality signs where assignments should be used.

Example:  You have x=x+1; That surely should be x:=x+1;
There are many such errors.

But there are other errors.

How did you do this by hand? By separation of variables? How?

@THAPELO You need to give the parameters to dsolve in one way or the other:

1. By assigning to the parameters before using dsolve, as in b:=1;  etc.
Be aware though that 'gamma' is Euler's constant in Maple, so you won't be allowed to assign to gamma. You could just call it something else, like gamma1.

2. A better way (in my opinion) is to define a set of equations (not assignments):

params:={b=1, d=0.2, k=1, alpha=0.5, gamma=0.25, mu=0.3};

and then use   eval({sys1, sys2, sys3},params)  as your concrete system. This is what I do below.
Using gamma is OK here as you are not assigning to it.

3. If you experiment a lot with varying parameters you may want to use the option parameters in dsolve

dsolve( {...ivp...}, numeric, parameters=[b, d, k, alpha, gamma, mu]);

This needs some explaining, which I shall not do now, but can give you an example if you wish.

##
If some parameters are fixed during the whole session you could assign to those as in point 1 and then use method 2 or 3 for the rest.

##############
restart;
sys1:= diff(S1(t),t) = b - d*S1(t) - k*S1(t)*S2(t)/(1+ alpha*S1(t)) - gamma*S3(t);
sys2:= diff(S2(t),t) = k*S1(t)*S2(t)/(1+ alpha*S1(t)) - (d + mu)*S2(t);
sys3:= diff(S3(t),t) = mu*S2(t) - (d + gamma)*S3(t);
##Method 2:
params:={b=1, d=0.2, k=1, alpha=0.5, gamma=0.25, mu=0.3};

Sol1:=dsolve(eval({sys1, sys2, sys3},params) union {S1(0)=1.35, S2(0)=0.9, S3(0)=0.45}, numeric);

plots:-odeplot(Sol1,[[t,S1(t)],[t,S2(t)],[t,S3(t)]],0..25);
#Values at t=25:
Sol1(25);
##Equilibria:
solve(rhs~({sys1, sys2, sys3}),{S1(t),S2(t),S3(t)});
eval([%],params);
map(subs,%,[S1(t),S2(t),S3(t)]);

@mathsstudent93 Just one more observation:

If you impose boundary conditions for s and/or D[1](s) at both ends of the interval you seem to get a trustworthy solution also for your original setup.

@mathsstudent93 I don't have any idea about what is going on, but do have a couple of observations:

1. You are missing a multiplication sign in DS just after (C2+3000.0).
  This is not the problem though.

2. There are some large and small numbers in the problem. So I tried a very simple version where every coefficient appearing is 1.
The code used was:

restart;
DS := diff(s(x, t), t) = diff(s(x, t), x, x)+b(x, t)/(1+b(x, t))-s(x, t)-e(x, t);
DE := diff(e(x, t), t) = s(x, t)-e(x, t)-e(x, t)*r(x, t)+b(x, t);
DR := diff(r(x, t), t) = 1+b(x, t)/(1+b(x, t))-r(x, t)-e(x, t)*r(x, t)+b(x, t);
DB := diff(b(x, t), t) = -b(x, t)+e(x, t)*r(x, t);
sys := DS, DE, DR, DB:
ivp := e(x, 0) = 1, s(x, 0) = 1, r(x, 0) = 1, b(x, 0) = 1, D[1](s)(0, t) = 0, s(0, t) = 1;
pde := pdsolve({sys}, {ivp}, numeric, time = t, spacestep=0.1, range = 0 .. 1);
pde:-plot(b(x, t), t = 1);
pde:-animate(b,t=1);
pde:-animate(s,t=1);



#OK you do get something out of there, but try changing spacestep to 0.01.

I suggest giving us the equations (as text) including (of course) initial and boundary conditions.

Notice that the help page for pdsolve/numeric says

 "The first mode of operation uses the default method, which is a centered implicit scheme, and is capable of finding solutions for higher order PDE or PDE systems. The PDE systems must be sufficiently close to a standard form for the method to find the numerical solution.
- The second mode of operation is strictly educational, and allows specification of a particular method to solve a PDE. This mode is restricted to a single PDE."

Since you have a system of 4 pdes you don't have any choice.

@J4James I suppose you want Maple to solve the last example.

Here I assume that past history is y(z)=1 for z<=0. (I don't know a way of giving dsolve a non-constant history as things are now with an alpha-version).

restart;
dde:=diff(y(z),z)=1-a*y(z)*y(z-xi)^n/(1+y(z-xi)^n);
params:={a=5,xi=0.6,n=4};
res:=dsolve({eval(dde,params),y(0)=1},numeric,range=0..100);
plots:-odeplot(res,[z,y(z)],0..100,refine=1,size=[900,200]);

You can remove the range option from dsolve and consequently the refine option from odeplot if you like.

@ribragimov Sysq has to be what it is in your worksheet and what I have above because then everything but delta1 and u1 are known.
The reason for finding the second derivatives of u0 and delta0 (via EQS) is that in order to solve the system Sysq and Sysqq we need to differentiate Sysq since it is not a differential equation. This can be done either explicitly (as I do) or by dsolve which will use a DAE-method (see the help under dsolve,numeric).
Since the system doesn't know how to differentiate F,F1,Fdiff,F1diff it will complain, however.
Rather than defining differentiation rules for those functions, I chose to do the differentiation of Sysq myself. That is how Sysqdiff came about. Notice that it only has u1 and delta1 as unknowns. 

With Sysq and Sysqdiff dsolve won't have to use a DAE-method, but will just use the default rkf45-method.

@Lal You can use piecewise:

f:=x->piecewise(x<Pi,exp(x),x^(-2));
plot(f(x),x=0..2*Pi,discont=true);
a:=0: b:=2*Pi: p:=b-a:
fp:=f(x-floor( (x-a)/p)*p);
plot(fp,x=-4*Pi..4*Pi,discont=true);