I tried your two equations with b = 1.

I got the same result in the two cases.

I believe it is because Maple starts by solving for f''', which involves squaring both sides.

@mah00 The system is not autonomous, i.e. the right hand sides depend explicitly on t (not only on x(t) and y(t)).

In that case arrows cannot be drawn since they would not depend on (x, y) only.

The wrapping in 'evalc' makes f1 and g1 into explicit real valued functions of the real variable t before being given to DEplot. Apparently DEplot doesn't handle Re and Im very well, if at all. 

Try e.g. the simple example

DEplot(diff(x(t),t)=evalc(Re(-x(t)+I*sin(t))),x(t),t=0..5,[x(0)=1]);

and try removing 'evalc'. You may want to look at the help for 'evalc'.

@mah00 The system is not autonomous, i.e. the right hand sides depend explicitly on t (not only on x(t) and y(t)).

In that case arrows cannot be drawn since they would not depend on (x, y) only.

The wrapping in 'evalc' makes f1 and g1 into explicit real valued functions of the real variable t before being given to DEplot. Apparently DEplot doesn't handle Re and Im very well, if at all. 

Try e.g. the simple example

DEplot(diff(x(t),t)=evalc(Re(-x(t)+I*sin(t))),x(t),t=0..5,[x(0)=1]);

and try removing 'evalc'. You may want to look at the help for 'evalc'.

I had more success with abserr=0.1e-5 (the other choices unaltered).

I also tried the default method with maxfun=0:
dsolve(dsys union isc, numeric, maxfun=0);
There wasn't much difference in the graphs:

plots:-odeplot(dsol,[[t,P(t)],[t,S(t)],[t,ES(t)],[t,E(t)]],0..100,thickness=[1,2,3,4]);

I had more success with abserr=0.1e-5 (the other choices unaltered).

I also tried the default method with maxfun=0:
dsolve(dsys union isc, numeric, maxfun=0);
There wasn't much difference in the graphs:

plots:-odeplot(dsol,[[t,P(t)],[t,S(t)],[t,ES(t)],[t,E(t)]],0..100,thickness=[1,2,3,4]);

@Markiyan Hirnyk You may consult

http://en.wikipedia.org/wiki/Weak_solution

which happens to have the pde in question as an example. It is easily verified that the solution I gave above satisfies (2) in the Wikipedia article.

@Markiyan Hirnyk You may consult

http://en.wikipedia.org/wiki/Weak_solution

which happens to have the pde in question as an example. It is easily verified that the solution I gave above satisfies (2) in the Wikipedia article.

@Markiyan Hirnyk My point with this example is that the numerical algorithm used by Maple provides a solution consistent with the one found by the exact solver after applying the initial and boundary conditions.

Clearly that solution is not continuous, so the question is really: What kind of solutions are you allowing? Should the solution solve the pde for all (z,t) in [0,1]x[0,1] or only for almost every (z,t) in [0,1]x[0,1]?

@Markiyan Hirnyk My point with this example is that the numerical algorithm used by Maple provides a solution consistent with the one found by the exact solver after applying the initial and boundary conditions.

Clearly that solution is not continuous, so the question is really: What kind of solutions are you allowing? Should the solution solve the pde for all (z,t) in [0,1]x[0,1] or only for almost every (z,t) in [0,1]x[0,1]?

The initial conditions are allowed to be discontinuous at a point:

eq:=diff(pA(z,t),t)=-diff(pA(z,t),z);

ibc:= [pA(z,0) = 0, pA(0,t) = 1];

res:=pdsolve(eq,ibc,numeric,time=t,range=0..1,spacestep=1e-3);

res:-animate(z=0..1,t=0..1,numpoints=300);
res:-plot3d(z=0..1,t=0..1,axes=boxed);

This simple equation is easily solved:

pdsolve(eq);
with result pA(z, t) = _F1(-z+t). The conditions ibc then imply

pA(z,t) = piecewise(t<z,0,1) (i.e. Heaviside(t-z))

plot3d(piecewise(t<z,0,1),z=0..1,t=0..1,axes=boxed);

The initial conditions are allowed to be discontinuous at a point:

eq:=diff(pA(z,t),t)=-diff(pA(z,t),z);

ibc:= [pA(z,0) = 0, pA(0,t) = 1];

res:=pdsolve(eq,ibc,numeric,time=t,range=0..1,spacestep=1e-3);

res:-animate(z=0..1,t=0..1,numpoints=300);
res:-plot3d(z=0..1,t=0..1,axes=boxed);

This simple equation is easily solved:

pdsolve(eq);
with result pA(z, t) = _F1(-z+t). The conditions ibc then imply

pA(z,t) = piecewise(t<z,0,1) (i.e. Heaviside(t-z))

plot3d(piecewise(t<z,0,1),z=0..1,t=0..1,axes=boxed);

'dsolve' should be 'pdsolve' as in the original question. I must have copied it wrong. I have corrected it above.

'dsolve' should be 'pdsolve' as in the original question. I must have copied it wrong. I have corrected it above.

The definition of B does not agree with ode[1], and in the definition of srPP you probably want srP instead of hrP.

Do you have any reason whatsoever why s1 and s2 (as defined last!) has anything to do with the difference between the unknown exact solution and the one found numerically by the method rk45 with abserr=10^(-10),relerr=10^(-10)??