Although this method is very fast, it seems to run into problems when the degree is 10000:

restart;
F:=randpoly(epsilon,degree=10000,dense):
simplify(F, {epsilon^2=1});

Error, (in solve_eqn) too many levels of recursion

Added: It seems that the revised version SE2 of my proposed procedure SE above is faster than using simplify/siderels.

@koonduck Yes, you can use odeplot like this

odeplot(dsol,[x,theta(x)],0..40);

# or you can use plot since in your worksheet you pulled out the theta-procedure in the line
dsol1 := subs(dsol, theta(x));
plot(dsol1,0..40);

@koonduck Yes, you can use odeplot like this

odeplot(dsol,[x,theta(x)],0..40);

# or you can use plot since in your worksheet you pulled out the theta-procedure in the line
dsol1 := subs(dsol, theta(x));
plot(dsol1,0..40);

You may try using a smaller number of steps (or a smaller stepsize), then you will see wiggling.

Your system is similar to a model for the damped pendulum (except for the factor 1/1.03 in the sine function).

Thus I would have expected to see damping. However, I don't.

Now try Maple's own Euler method for comparison with your version and also with Maple's default numerical method rk45:

eq:=diff(x(t),t,t)=f2(t,x(t),diff(x(t),t));
resEuler:=dsolve({eq,x(0)=5,D(x)(0)=0},numeric,method=classical[foreuler],stepsize=.05);
plots:-odeplot(resEuler,[t,x(t)],0..5);
#Now the default Runge-Kutta-Fehlberg method:
resRKF45:=dsolve({eq,x(0)=5,D(x)(0)=0},numeric,stepsize=.05);
plots:-odeplot(resRKF45,[t,x(t)],0..5);

I guess this says something about Euler's method!

You may try using a smaller number of steps (or a smaller stepsize), then you will see wiggling.

Your system is similar to a model for the damped pendulum (except for the factor 1/1.03 in the sine function).

Thus I would have expected to see damping. However, I don't.

Now try Maple's own Euler method for comparison with your version and also with Maple's default numerical method rk45:

eq:=diff(x(t),t,t)=f2(t,x(t),diff(x(t),t));
resEuler:=dsolve({eq,x(0)=5,D(x)(0)=0},numeric,method=classical[foreuler],stepsize=.05);
plots:-odeplot(resEuler,[t,x(t)],0..5);
#Now the default Runge-Kutta-Fehlberg method:
resRKF45:=dsolve({eq,x(0)=5,D(x)(0)=0},numeric,stepsize=.05);
plots:-odeplot(resRKF45,[t,x(t)],0..5);

I guess this says something about Euler's method!

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.

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