Now what I did is:

E := (1/2)*(1/rho1^2+1/rho1^2)-1/sqrt(z1^2+rho1^2)-1/sqrt((R-z1)^2+rho1^2)-1/sqrt(z1^2+rho1^2)-1/sqrt((R-z1)^2+rho1^2)+1/R+1/sqrt((2*rho1)^2+(R-2*z1)^2)

Er1 := diff(E, rho1) = 0

Ez1 := diff(E, z1) = 0

sol := proc (r) options operator, arrow; fsolve(eval({Er1, Ez1}, R = r), {z1 = -1 .. 3, rho1 = 0 .. 3}) end proc

sol(.5)

solE := proc (s) options operator, arrow; eval(E, {op(1, sol(s)), op(2, sol(s)), R = s}) end proc

solE(.5)

plot(proc (s) options operator, arrow; solE(s) end proc, 0 .. 5, view = -.5 .. -1.2)

For me at 2.3/2.4 the plot is not continuous.

I think this is due to the fsolve command mixing up different solutions.

There is one solution which behaves like shown there in the plot for large s but also which continues less steep

e.g. if I specify rho1=0..1.02 you see the continuous function.

But how can I select the solution whose behaviour for large s you see when increasing the search bound. e.g. rho1=0..3

@Preben Alsholm 

You wrote:

Actually, the indicial equation returned by 
indicialeq(ode2, r, 0, X(y));
r^3-3*r^2+2*r = 0
is wrong!
The indicial equation is
-10+(33/2)*r-(15/2)*r^2+r^3 = 0;

From me? I always had the second equation, so the roots 1,5/2,4

how did you get the first equation?

@Preben Alsholm 

You wrote:

Actually, the indicial equation returned by 
indicialeq(ode2, r, 0, X(y));
r^3-3*r^2+2*r = 0
is wrong!
The indicial equation is
-10+(33/2)*r-(15/2)*r^2+r^3 = 0;

From me? I always had the second equation, so the roots 1,5/2,4

how did you get the first equation?

@Preben Alsholm I tried a little bit. with the right a4 and a5 initial condition 

dsolve(ode2, X(y), 'series', y = 0)

is the same as your truncated solutions sol2 and sol3 respectively

But I still dont understand what sort of solution 

dsolve(ode2, X(y), 'formal_series', y = 0); gives. I know its the expansion about the corresponding singularity. But I dont get why its a polynomial. I thought the other 2 solutions were no polynomials. And it doesnt say that it ignores higher orders like usual when doing series stuff.

You are talking about that y=0 is ignore since it only considers 'coeffs'=whatever only.

if I take 'coeffs'=polynomial than it looks for a solution which is polynomial in the coefficients or what?

same for rational?

In the end:

When having a recurrence like this  (which at first is valid for -infinity < n < infinity, right?) you have to look which parts of the recurrence seperate from the other. like:n=-5..4

520 z a(-5) + (-384 z - 384 + 12 lambda) a(-4) + 308 a(-3) = 0
308 z a(-4) + (-210 z - 210 + 10 lambda) a(-3) + 162 a(-2) = 0
162 z a(-3) + (-96 z - 96 + 8 lambda) a(-2) + 70 a(-1) = 0
70 z a(-2) + (-30 z - 30 + 6 lambda) a(-1) + 20 a(0) = 0
20 z a(-1) + 4 lambda a(0) = 0
(6 z + 2 lambda + 6) a(1) - 2 a(2) = 0
-2 a(1) z + 2 a(3) = 0
2 a(2) z + (-6 z - 6 - 2 lambda) a(3) = 0
-4 lambda a(4) - 20 a(5) = 0
-20 z a(4) + (30 z + 30 - 6 lambda) a(5) - 70 a(6) = 0

there is no relation between a(3) and a(4) so it truncates for n>3

What about n<0 possibly I could have a solution from -infinity..3, so a laurent expansion? or am I on the wrong way :-(

@Preben Alsholm I tried a little bit. with the right a4 and a5 initial condition 

dsolve(ode2, X(y), 'series', y = 0)

is the same as your truncated solutions sol2 and sol3 respectively

But I still dont understand what sort of solution 

dsolve(ode2, X(y), 'formal_series', y = 0); gives. I know its the expansion about the corresponding singularity. But I dont get why its a polynomial. I thought the other 2 solutions were no polynomials. And it doesnt say that it ignores higher orders like usual when doing series stuff.

You are talking about that y=0 is ignore since it only considers 'coeffs'=whatever only.

if I take 'coeffs'=polynomial than it looks for a solution which is polynomial in the coefficients or what?

same for rational?

In the end:

When having a recurrence like this  (which at first is valid for -infinity < n < infinity, right?) you have to look which parts of the recurrence seperate from the other. like:n=-5..4

520 z a(-5) + (-384 z - 384 + 12 lambda) a(-4) + 308 a(-3) = 0
308 z a(-4) + (-210 z - 210 + 10 lambda) a(-3) + 162 a(-2) = 0
162 z a(-3) + (-96 z - 96 + 8 lambda) a(-2) + 70 a(-1) = 0
70 z a(-2) + (-30 z - 30 + 6 lambda) a(-1) + 20 a(0) = 0
20 z a(-1) + 4 lambda a(0) = 0
(6 z + 2 lambda + 6) a(1) - 2 a(2) = 0
-2 a(1) z + 2 a(3) = 0
2 a(2) z + (-6 z - 6 - 2 lambda) a(3) = 0
-4 lambda a(4) - 20 a(5) = 0
-20 z a(4) + (30 z + 30 - 6 lambda) a(5) - 70 a(6) = 0

there is no relation between a(3) and a(4) so it truncates for n>3

What about n<0 possibly I could have a solution from -infinity..3, so a laurent expansion? or am I on the wrong way :-(

@digerdiga Why do I get when using:

dsolve(ode, X(y), 'formal_series', y = 0)

3 polynomial solutions 

but when using:

dsolve(ode, X(y), 'series', y = 0)

None of the solutions  truncate..

Thats weird

@digerdiga Why do I get when using:

dsolve(ode, X(y), 'formal_series', y = 0)

3 polynomial solutions 

but when using:

dsolve(ode, X(y), 'series', y = 0)

None of the solutions  truncate..

Thats weird

@Preben Alsholm Thanks. Thats what I wanted to know ;)

How can I check when calculating the recurrence what specific indicial root it corresponds to?

(other than knowing that the solution is of the form y*(a0+a1*y+....) so it would be 1)

Ah and why does the recurrence say division by zero at some value?

When having a polynomial wouldnt the recurrence truncate at this specific value i.e. would be zero for the following values...and not "divergent"

@Preben Alsholm Thanks. Thats what I wanted to know ;)

How can I check when calculating the recurrence what specific indicial root it corresponds to?

(other than knowing that the solution is of the form y*(a0+a1*y+....) so it would be 1)

Ah and why does the recurrence say division by zero at some value?

When having a polynomial wouldnt the recurrence truncate at this specific value i.e. would be zero for the following values...and not "divergent"

@Preben Alsholm does it mean that this solution does not exist when it says division by zero?

ode := (25*x^2-15*x*z+x*lambda-15*x+10*z)*(diff(X(x), x))+(-9*x^2+(15/2)*x*z+(21/2)*x^3-9*x^2*z)*(diff(X(x), x, x))+(x^4-x^3*z-x^3+x^2*z)*(diff(X(x), x, x, x))+10*X(x)*x+2*X(x)*lambda = 0;

with(PDEtools);

ode2 := collect(radnormal(y*dchange(x = 1/y, ode)), diff);

with(DEtools);

indicialeq(ode2, y, 0, X(y));

solve(%);

with(gfun);

f := diffeqtorec(ode2, X(y), a(n)); subs(n = m-2, f);

a := proc (m) -((-2*z*(m-2)^3+9*z*(m-2)^2-9*z*(m-2))*a(m-2)+((2*z+2)*(m-2)^3+(-6*z-6)*(m-2)^2+(-2*z-2*lambda-2)*(m-2)+6*z+2*lambda+6)*a(m-1))/(-2*(m-2)^3+3*(m-2)^2+3*m-8) end proc;

a(-1) := 0; a(0) := 1

a(1)

Error, (in a) numeric exception: division by zero

@Preben Alsholm does it mean that this solution does not exist when it says division by zero?

ode := (25*x^2-15*x*z+x*lambda-15*x+10*z)*(diff(X(x), x))+(-9*x^2+(15/2)*x*z+(21/2)*x^3-9*x^2*z)*(diff(X(x), x, x))+(x^4-x^3*z-x^3+x^2*z)*(diff(X(x), x, x, x))+10*X(x)*x+2*X(x)*lambda = 0;

with(PDEtools);

ode2 := collect(radnormal(y*dchange(x = 1/y, ode)), diff);

with(DEtools);

indicialeq(ode2, y, 0, X(y));

solve(%);

with(gfun);

f := diffeqtorec(ode2, X(y), a(n)); subs(n = m-2, f);

a := proc (m) -((-2*z*(m-2)^3+9*z*(m-2)^2-9*z*(m-2))*a(m-2)+((2*z+2)*(m-2)^3+(-6*z-6)*(m-2)^2+(-2*z-2*lambda-2)*(m-2)+6*z+2*lambda+6)*a(m-1))/(-2*(m-2)^3+3*(m-2)^2+3*m-8) end proc;

a(-1) := 0; a(0) := 1

a(1)

Error, (in a) numeric exception: division by zero

@digerdiga btw isn't it possible to obtain the 2 other roots of the indicial equation?. i.e. the expansion is of the form

X(x):=sum_{n=0}^{infinity} a_n x^{n+r}

where r is such that after shifting the sums x^{r-1} and x^{r-2} vanish

for each root there is a solution so there should be 3 recurrence relations in my situation.

How do I know which root diffeqtorec chose?

@digerdiga btw isn't it possible to obtain the 2 other roots of the indicial equation?. i.e. the expansion is of the form

X(x):=sum_{n=0}^{infinity} a_n x^{n+r}

where r is such that after shifting the sums x^{r-1} and x^{r-2} vanish

for each root there is a solution so there should be 3 recurrence relations in my situation.

How do I know which root diffeqtorec chose?

@Markiyan Hirnyk well it seems to be due to the specific worksheet I think.

In another worksheet it does work.

I just used:

with(gfun);
diffeqtorec(odeX1, X(x), a(n));
Error, (in gfun:-diffeqtorec) invalid input: getname expects its 1st argument, yofz, to be of type function(name), but received -(3/2)/x(n)+(3/2)/(x(n)-1)+(3/2)/(x(n)-z(n))

odeX1 is the name I called it.

@Markiyan Hirnyk well it seems to be due to the specific worksheet I think.

In another worksheet it does work.

I just used:

with(gfun);
diffeqtorec(odeX1, X(x), a(n));
Error, (in gfun:-diffeqtorec) invalid input: getname expects its 1st argument, yofz, to be of type function(name), but received -(3/2)/x(n)+(3/2)/(x(n)-1)+(3/2)/(x(n)-z(n))

odeX1 is the name I called it.