@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);

@Lal p is just the length of the interval a..b. The code is written so it can be used for any function h defined on a..b.

The argument  x-floor( (x-a)/p)*p) will lie in the interval a..b for any real x. To be precise it will lie in the half-open interval [a,b).
Try
plot(x-floor( (x-a)/p)*p,x=-4*Pi..4*Pi,scaling=constrained,discont=true);

You may also want to look at the help page for floor:
?floor

@sssMMM To continue where I left off after having defined Fdiff and F1diff I did:

Sysq := ((9/10)*F1(x)*F(x)*F1diff(x)-(3/10)*Fdiff(x)*F1(x)*F1(x)+(3/2)*F1(x))*F(x) = (1/5)*F(x)*delta1(x)-(3/10)*u1(x);
Sysqq := u1(x)*Fdiff(x)+F1(x)*(diff(delta1(x), x))+delta1(x)*F1diff(x)+F(x)*(diff(u1(x), x)) = 0;

sbs1:={diff(u0(x),x)=F1diff(x),diff(delta0(x),x)=Fdiff(x),delta0(x)=F(x),u0(x)=F1(x)};
Fdiff2:=unapply(subs(diff(EQS,x),sbs1,diff(delta0(x),x,x)),x,numeric);
F1diff2:=unapply(subs(diff(EQS,x),sbs1,diff(u0(x),x,x)),x,numeric);

sbs2:={diff(F(x),x)=Fdiff(x),diff(F1(x),x)=F1diff(x),diff(Fdiff(x),x)=Fdiff2(x),diff(F1diff(x),x)=F1diff2(x)};
Sysqdiff:=subs(sbs2,diff(Sysq,x));

sol:=dsolve({Sysqdiff,Sysqq,delta1(0) = 10^(-8), u1(0) = 10^(-8)},numeric,known=[F,F1,Fdiff,F1diff,Fdiff2,F1diff2]);
odeplot(sol,[[x,u1(x)],[x,delta1(x)]],0..10);

@J4James I think that Allan Wittkopf's second reply answers that completely.

What is in that html file you have at the bottom? I cannot access it.

Try again, but this time after a restart.
Thus do

restart;
int(x^2+3*x, x);

@J4James Your headline is "Method of Steps" and the content is "I'm wondering maybe Maple can still handle DDEs with a Proc.".

I fail to understand your question.

############
Do you mean:
How does the procedure (when using procedural input) have to look for a delay differential equation if it is at all available to the user?
As you will understand from Allan Wittkopf's reply the feature is still under development and is so far undocumented.
So I don't know how it is or how it will become.

To make clear what I'm trying to say. Consider a usual ode:
ode:=diff(x(t),t,t)+t*x(t)=0;
sol:=dsolve({ode,x(0)=1,D(x)(0)=0},numeric);
plots:-odeplot(sol,[t,x(t)],0..15);
#Now use procedural input:
p:=proc(n,t,Y,YP)
  YP[1]:=Y[2];
  YP[2]:=-t*Y[1]
end proc;
res:=dsolve(numeric,number=2,procedure=p,start=0,initial=Array([1,0]),procvars=[x(t),diff(x(t),t)]);
plots:-odeplot(res,[t,x(t)],0..15);

So is there an analogue or will there be?
There is no point in me guessing.

@Allan Wittkopf Thanks for the info. Good to hear!

@J4James My version is Maple 18.02.

And sorry for the confusing edits, but it came as a real surprise to me that Maple 18 can handle DDEs.

I shall experiment with this, but I just noticed that it appears to handle time-dependent delays too.

@J4James I didn't see the problem, since I just let Maple 18 do it it as stated.
It replied with [Length of output exceeds limit of 1000000], which didn't disturb me. I proceeded with
plots:-odeplot(res,[t,x(t)],0..100);

which produced the following

This morning somewhat more clear-headed (apparently not much more) I noticed the delay, and as I was pretty sure that dsolve cannot handle DDEs I was sure something was wrong.
I tried your system on my own ddesolve assuming that the history for t=-10..0 is given by the constant values in the ics (tell me if the history should be something else).
The graph produced for x on the same interval was to the eye exactly the same!

Then I tried Maple 17 and I got the same error as reported by you.

Conclusion:
This feature is new in Maple 18.

As far as I know this is not documented. If I'm wrong I should like to know.

Try this simple example:
restart;
ode:=diff(x(t), t)=x(t-1);
infolevel[`dsolve/numeric`]:=1;
res:=dsolve({ode,x(0)=1},numeric);
plots:-odeplot(res,[t,x(t)],t=0..5);
#The assumption made by dsolve/numeric apparently is that the history is given by x(t)=1 for -1<t<=0.

Could you show us what you did with text, please, not images.
I'm confused about your use of x,y,x1,y1 in your vector equation. Which represents x' and y'?