@escorpsy It is difficult to make nice looking textbook like pictures.
But here is a rough worksheet which initially uses DEplot3d and at the end odeplot.
Explanations could be much more plentiful, but that takes time.

ODEsys12-12-14.mw

@escorpsy It is difficult to make nice looking textbook like pictures.
But here is a rough worksheet which initially uses DEplot3d and at the end odeplot.
Explanations could be much more plentiful, but that takes time.

ODEsys12-12-14.mw

Maybe you could give us an example? Because the following works:

p:=randpoly(x);
q:=eval(p,x=sqrt(y))+a*sqrt(y);
collect(q,y^(1/2));

Added: OK, Patrick got there before I did, but I'll leave it.

@Thomas Richard I added a convert/polynom to your lines. That clearly exhibits the problem as a typesetting bug.

restart;
interface(typesetting);
interface(typesetting=extended): # bad
s1:=series(1/epsilon-log(4),epsilon=0);
whattype(s1);
convert(s1,polynom);
interface(typesetting=standard): # good
s2:=series(1/epsilon-log(4),epsilon=0);
whattype(s2);
convert(s2,polynom);

@J4James Yes, there are negative eigenvalues, but for another of the solutions for f. The boundary value problem for f has (at least) two solutions.

I have revised yet again the uploaded worksheet.

@J4James Yes, there are negative eigenvalues, but for another of the solutions for f. The boundary value problem for f has (at least) two solutions.

I have revised yet again the uploaded worksheet.

@J4James Do you have any evidence of negative eigenvalues? Is the system considered (and the constants) in the image you gave exactly the same as yours? I noticed a difference for the image of the positive eigenvalues.

I have not been able to find any evidence myself although I have tried. I have uploaded an enlarged version of the worksheet, which also corrects the error statement in the procedure L ('tol' corrected to 'tolerance').

@J4James Do you have any evidence of negative eigenvalues? Is the system considered (and the constants) in the image you gave exactly the same as yours? I noticed a difference for the image of the positive eigenvalues.

I have not been able to find any evidence myself although I have tried. I have uploaded an enlarged version of the worksheet, which also corrects the error statement in the procedure L ('tol' corrected to 'tolerance').

@PatrickT What did you intend the last line to do?

Try

seq( Q(b,Lambda), params = [b=1,Lambda=1] ) ;

@PatrickT What did you intend the last line to do?

Try

seq( Q(b,Lambda), params = [b=1,Lambda=1] ) ;

@escorpsy Don't you just mean the projection of the orbits onto the xy-plane?

If so, you can do as indicated in an earlier response, i.e. use DEplot (not DEplot3d), but add the scene option: scene=[x,y]. Beware though that there likely will be problems because some orbits fly off to infinity thereby creating lots of data with HFloat(undefined). It helps to make the timerange shorter in both ends and make the stepsize smaller (0.001).
Added: Because of the absence of similar problems with DEplot3d it may be a better idea simply to rotate the 3d plot produced by DEplot3d. This can also be done programmatically by using the orientation option. Suppose the plot generated by DEplot3d is assigned to the variable p, then
display(p,bg,orientation=[90,0,180]); #xy-plane
display(p,bg,orientation=[90,90,180]);#xz-plane
display(p,bg,orientation=[90,90,90]);##zy-plane



Another approach is to use odeplot instead as indicated in the reference given by Patrick. I haven't looked into that, but with the same definitions as before you can try:

ptsxy:=subsop~(3=NULL,pts);
bgxy:=pointplot(ptsxy,symbol=solidcircle,symbolsize=20,color=blue):
RES:=dsolve({sys,x(0)=x00,y(0)=y00,z(0)=z00},numeric,parameters=[x00,y00,z00]);
pp:=proc(x0,y0,z0,timerange)
  RES(parameters=[x0,y0,z0]);
  plots:-odeplot(RES,[x(t),y(t)],timerange, _rest)
end proc:
# pp accepts all optional arguments to odeplot (since they are passed directly to odeplot via '_rest').
pp(.4,2,1,-0.9..0,view=[-5..5,-5..5]);
PL:=[seq(seq(seq(pp(x0+i*dx/N,y0+j*dy/N,z0+k*dz/N,-1..1,view=[x0..x0+dx,y0..y0+dy]),k=1..N-1),j=1..N-1),i=1..N-1)]:
display(bgxy,PL);

@escorpsy Don't you just mean the projection of the orbits onto the xy-plane?

If so, you can do as indicated in an earlier response, i.e. use DEplot (not DEplot3d), but add the scene option: scene=[x,y]. Beware though that there likely will be problems because some orbits fly off to infinity thereby creating lots of data with HFloat(undefined). It helps to make the timerange shorter in both ends and make the stepsize smaller (0.001).
Added: Because of the absence of similar problems with DEplot3d it may be a better idea simply to rotate the 3d plot produced by DEplot3d. This can also be done programmatically by using the orientation option. Suppose the plot generated by DEplot3d is assigned to the variable p, then
display(p,bg,orientation=[90,0,180]); #xy-plane
display(p,bg,orientation=[90,90,180]);#xz-plane
display(p,bg,orientation=[90,90,90]);##zy-plane



Another approach is to use odeplot instead as indicated in the reference given by Patrick. I haven't looked into that, but with the same definitions as before you can try:

ptsxy:=subsop~(3=NULL,pts);
bgxy:=pointplot(ptsxy,symbol=solidcircle,symbolsize=20,color=blue):
RES:=dsolve({sys,x(0)=x00,y(0)=y00,z(0)=z00},numeric,parameters=[x00,y00,z00]);
pp:=proc(x0,y0,z0,timerange)
  RES(parameters=[x0,y0,z0]);
  plots:-odeplot(RES,[x(t),y(t)],timerange, _rest)
end proc:
# pp accepts all optional arguments to odeplot (since they are passed directly to odeplot via '_rest').
pp(.4,2,1,-0.9..0,view=[-5..5,-5..5]);
PL:=[seq(seq(seq(pp(x0+i*dx/N,y0+j*dy/N,z0+k*dz/N,-1..1,view=[x0..x0+dx,y0..y0+dy]),k=1..N-1),j=1..N-1),i=1..N-1)]:
display(bgxy,PL);

@wdarrel Applying the Chain rule (not Maple syntax):

diff(u,t) = diff(w,tau)*diff(tau,t) + diff(w,y)*diff(y,t) = diff(w,tau)*1/2 + diff(w,y)*(-1/2)

diff(u,x) = diff(w,tau)*diff(tau,x) + diff(w,y)*diff(y,x) = diff(w,tau)*1/2 + diff(w,y)*1/2

Adding those you see the cancellation of the diff(w,y) terms.

In Maple and first in general:

restart;
#The variable change:
V:={y=f(x,t),tau=g(x,t)};
eq0:=u(x,t)=w(y,tau);
eq:=subs(V,eq0);
#Differentiating:
eqt:=diff(eq,t);
eqx:=diff(eq,x);
#Finding the variable change in our case, i.e. finding y = f(x,t) and tau=g(x,t):
S:=solve({ x=tau+y,t=tau-y },{tau,y});
f:=(x,t)->(x+t)/2;
g:=(x,t)->(x-t)/2;
eqx;
eqt;
%%+%;
subs({x=tau+y,t=tau-y},rhs(%));
convert(%,diff);

@wdarrel Applying the Chain rule (not Maple syntax):

diff(u,t) = diff(w,tau)*diff(tau,t) + diff(w,y)*diff(y,t) = diff(w,tau)*1/2 + diff(w,y)*(-1/2)

diff(u,x) = diff(w,tau)*diff(tau,x) + diff(w,y)*diff(y,x) = diff(w,tau)*1/2 + diff(w,y)*1/2

Adding those you see the cancellation of the diff(w,y) terms.

In Maple and first in general:

restart;
#The variable change:
V:={y=f(x,t),tau=g(x,t)};
eq0:=u(x,t)=w(y,tau);
eq:=subs(V,eq0);
#Differentiating:
eqt:=diff(eq,t);
eqx:=diff(eq,x);
#Finding the variable change in our case, i.e. finding y = f(x,t) and tau=g(x,t):
S:=solve({ x=tau+y,t=tau-y },{tau,y});
f:=(x,t)->(x+t)/2;
g:=(x,t)->(x-t)/2;
eqx;
eqt;
%%+%;
subs({x=tau+y,t=tau-y},rhs(%));
convert(%,diff);

@User7843 The line

evalf((eval(phi,t=0)-eval(phi,t=3600))/2/Pi);

gives the number of periods of the sin function. The result is approximately 2262 on the interval I was using.
Try removing evalf. That will show what it does.
I'm using the fact that phi is a monotone (decreasing) function of t.