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)??

@kinesimario How is Y (capital Y) defined?

@kinesimario How is Y (capital Y) defined?

@kinesimario As you write it there is nothing wrong once you have corrected ic2 as explained by Jarekk.

But why use a stiff solver?

deq2 := diff(x(t),t,t) = -x(t)/(x(t)^2+y(t)^2)^(3/2), diff(y(t), t,t) = -y(t)/(x(t)^2+y(t)^2)^(3/2);

ic2 := x(0) = 1, D(x)(0) = 0, y(0) = 0, D(y)(0) = 1;

rosen := dsolve({ic2, deq2}, type = numeric, range = 0 .. 2*Pi);

plots:-odeplot(rosen,[x(t),y(t)],0..2*Pi);

You may find some fun in trying the following:

rosen1 := dsolve({ic2, deq2}, numeric);
plots:-odeplot(rosen1,[x(t),y(t)],0..2*Pi,scaling=constrained);
?dsolve,interactive,numeric
rosen1(initial);
rosen1(initial=[1.5,0,0,1]);
plots:-odeplot(rosen1,[x(t),y(t)],0..33,scaling=constrained);
rosen1(initial=[2,0,0,1]);
plots:-odeplot(rosen1,[x(t),y(t)],-500..500,scaling=constrained);