@Preben Alsholm Thank you so much! I'd been stuck on this for three weeks!

@Preben Alsholm 

Thanks for your reply! I will have to stay with the default numeric solver then.

If it's not too cheeky, could you have a quick look at what I've got so far? It's a very ugly equation I'm afraid.

Ds := 0.057;
ks := 0.005;

kr := 0.03;

ke := 0.005;
ka := 1.8*10^(-10);
kd := 0.12;
ki :=0.043;
dl := 0.4;
C2 := 8000.0;
E1 := (kd+ki)/ka;
S1 := ((dl+ke)*E1+3000.0*ki)/dl;
DS := diff(s(x, t), t) = Ds*(diff(s(x, t), x, x))+(C2+3000.0)(3020.0+1510.0/dl)*b(x, t)/(C2+b(x, t))-dl*(s(x, t)-e(x, t))-ks*s(x, t);

DE := diff(e(x, t), t) = dl*(s(x, t)-e(x, t))-ke*e(x, t)-ka*e(x, t)*r(x, t)+kd*b(x, t);

DR := diff(r(x, t), t) = 90.0+27300.0*b(x, t)/(6.32*10^5+b(x, t))-kr*r(x, t)-ka*e(x, t)*r(x, t)+kd*b(x, t);

DB := diff(b(x, t), t) = -ki*b(x, t)+ka*e(x, t)*r(x, t)-kd*b(x, t);
sys := DS, DE, DR, DB;

ivp := e(x, 0.) = E1, s(x, 0.) = S1, r(x, 0.) = 3000.0, b(x, 0.) = 3000.0, (D[1](s))(0., t) = 0., s(0., t) = 10*S1;
pde := pdsolve({sys}, {ivp}, numeric, time = t, range = 0. .. 20.0);

Which it looks like it's solved, but then if I try to plot anything I receive:

plotb1 := pde:-plot(b(x, t), t = 0.0);

Error, (in pdsolve/numeric/plot) unable to compute solution for t<HFloat(0.0): matrix is singular

which I was told may mean that the Newton Method wasn't converging, though not by too reliable a source.

Do you have any idea what's going on?