But the variable names have nothing to do with it though.

restart;
with(RandomTools):
s := evalf(Generate(list(rational(denominator = 30), 10)));
s1 := evalf(Generate(list(rational(denominator = 30), 10)));
restart;
with(RandomTools):
s1 := evalf(Generate(list(rational(denominator = 30), 10)));
s := evalf(Generate(list(rational(denominator = 30), 10)));

A version not using identify:

interface(displayprecision=5);
evalindets(k,numeric,evalf[18]):
evalindets(%,`^`,x->op(1,x)^convert(op(2,x),rational));

A version not using identify:

interface(displayprecision=5);
evalindets(k,numeric,evalf[18]):
evalindets(%,`^`,x->op(1,x)^convert(op(2,x),rational));

@Axel Vogt Maple agrees:

sum(2*k-1,k=1..n);
expand(%);

@Axel Vogt Maple agrees:

sum(2*k-1,k=1..n);
expand(%);

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.

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.