Did you try the code in my comment? It reproduces the graphs on the top of page 5 rather exactly. It uses equtions 13 and 14. By r^2 is meant the inner product of r with r, i.e. <x,y>.<x,y> = x^2 + y^2.

@kagestodder r = (x, y) :

V2:=unapply(V(x,y)-11/2000*(x^2+y^2),x,y);
p1:=plot(V(x,0),x=-15..15,-5..0,caption="Gravitational potential"):
p2:=plot(V2(x,0),x=-15..15,-5..0,caption="Pseudo-potential"):
plots:-display(Array([p1,p2]));

@kagestodder r = (x, y) :

V2:=unapply(V(x,y)-11/2000*(x^2+y^2),x,y);
p1:=plot(V(x,0),x=-15..15,-5..0,caption="Gravitational potential"):
p2:=plot(V2(x,0),x=-15..15,-5..0,caption="Pseudo-potential"):
plots:-display(Array([p1,p2]));

@williamov I confess to not knowing much about this. However, I tried to work backwards hoping to gain some insight:

restart;
Sum(HermiteH(n,x)*z^n/n!,n=0..infinity);
Heven:=convert(HermiteH(2*k,x),hypergeom) assuming k::nonnegint;
Hodd:=convert(HermiteH(2*k+1,x),hypergeom) assuming k::nonnegint;
G1:=Sum(Heven*z^(2*k)/(2*k)!,k=0..infinity)+Sum(Hodd*z^(2*k+1)/(2*k+1)!,k=0..infinity);
G2:=exp(2*x*z-z^2);
#Small test:
evalf(eval(G1,{z=1.2345,x=.567}));
evalf(eval(G2,{z=1.2345,x=.567}));

@williamov I confess to not knowing much about this. However, I tried to work backwards hoping to gain some insight:

restart;
Sum(HermiteH(n,x)*z^n/n!,n=0..infinity);
Heven:=convert(HermiteH(2*k,x),hypergeom) assuming k::nonnegint;
Hodd:=convert(HermiteH(2*k+1,x),hypergeom) assuming k::nonnegint;
G1:=Sum(Heven*z^(2*k)/(2*k)!,k=0..infinity)+Sum(Hodd*z^(2*k+1)/(2*k+1)!,k=0..infinity);
G2:=exp(2*x*z-z^2);
#Small test:
evalf(eval(G1,{z=1.2345,x=.567}));
evalf(eval(G2,{z=1.2345,x=.567}));

williamov wrote P(0) = 1, which changes the constants, but not the approach:

restart;
test4 := diff(P(x),x$2)-2*x*diff(P(x),x)+2*n*P(x)=0;
#test4 := {diff(P(x),x$2)-2*x*diff(P(x),x)+2*n*P(x), P(0) = 1, (D(P))(0) = 1};
res:=dsolve(test4);
series(rhs(res),x=0,2);
convert(%,polynom) assuming x>0;
solve(identity(%=1+x,x),{_C1,_C2});
eval(res,%);
convert(%,hypergeom);

williamov wrote P(0) = 1, which changes the constants, but not the approach:

restart;
test4 := diff(P(x),x$2)-2*x*diff(P(x),x)+2*n*P(x)=0;
#test4 := {diff(P(x),x$2)-2*x*diff(P(x),x)+2*n*P(x), P(0) = 1, (D(P))(0) = 1};
res:=dsolve(test4);
series(rhs(res),x=0,2);
convert(%,polynom) assuming x>0;
solve(identity(%=1+x,x),{_C1,_C2});
eval(res,%);
convert(%,hypergeom);

I suggested infinity as default so that the two argument version will act as before.  It should also be safe to overwrite the library version of applyrule with infinity as default.

I suggested infinity as default so that the two argument version will act as before.  It should also be safe to overwrite the library version of applyrule with infinity as default.

@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.