@joshkalon tomleslie gave you the code.

I supply it here too.

restart;
ode:=diff(delta(x),x,x)+(1/2*0.3*(1+x)^3-0.7)/(0.3*(1+x)^3+0.7)/(1+x)*diff(delta(x),x)-3/2*0.3*(1+x)/(0.3*(1+x)^3+0.7)*delta(x)=0;
ics:=delta(0)=1,D(delta)(0)=-0.5;
dsolve({ode,ics});
res:=value(%) assuming x>0; #Presumably x>0 ?
eval(res,x=19); #Value at x=19
evalf(%); #In decimals
#Numerical check using numerical solution of the ode:
resnum:=dsolve({ode,ics},numeric);
plots:-odeplot(resnum,[x,delta(x)],0..19); #Plot of solution up till x=19
resnum(19); #Compares favorably with the result obtained from the exact solution.

Even for a flat earth you have the problem:

restart;
with(plots):
surf:=3;
V := Vector([diff(x(t), t), diff(y(t), t)]);
N := Vector([0,1]);
V_reflect := -(1+CR)*LinearAlgebra:-DotProduct(V,N,conjugate=false)*N+V;
CR := .8;
deqs := diff(y(t), t, t) = -9.81, diff(x(t), t, t) = 0;
ics := (D(x))(0) = 2, (D(y))(0) = 0, x(0) = 0, y(0) = 3.5;
sol := dsolve({deqs, ics}, {x(t), y(t)}, numeric,
   events = [[y(t) = 3, [temp = diff(x(t), t), diff(x(t), t) = V_reflect[1], diff(y(t), t) = subs(diff(x(t), t) = temp, V_reflect[2])]]], output = listprocedure);
xanim,yanim:=op(subs(sol, [x(t),y(t)]));
p1 := plot(surf, x = 0 .. 10, color = black):
animate(pointplot, [[xanim(t), yanim(t)], symbol = solidcircle, symbolsize = 15, color = black], t = 0 .. 3, frames = 150,background=p1);

By just using the initial condition pdsolve comes up with an answer, which happens also to satisfy the boundary condtitions:

pdsolve([pde,bc[3]]);

comes up with u(x, t) = x.

@fanicia Your list uses a capital X (upper case). In acer's reply it says x=X.

@J4James In this particular example you could do it in this (not very elegant) way:

res:=combine(PDEtools:-dchange({Xsubs2,Ysubs2,psubs2},Eq,[X,Y,P]));
indets(res,specfunc(anything,exp));
ln~(%) assuming real;
`=`(op(%));
simplify(%/epsilon);

@Markiyan Hirnyk The code exactly as I gave works for me in Maple 18, 17, 16, and 15.
So which Maple version are you using?

I don't have Maple 2015 yet.

@z23 I just varied n. Try p(2), p(3) ,p(4), etc.

By only n-2 eigenvalues outside the x-interval I simply mean that for abs(x) >1 (the even case) there are n-2 real eigenvalues and consequently 2 imaginary ones.
Thus for n=50 (your case) I am saying that there are 48 real eigenvalues when abs(x) > 1 and 2 imaginary ones. For abs(x)<=1 all 50 eigenvalues are real.

@Axel Vogt I realize that my answers are rather colorless and could benefit from a certain amount of display of results. I shall try to keep your comment in mind when answering in the future.

@acer This gives me an opportunity to ask you acer why you consistently (it seems) use restart: and not restart;

Both work, yes, both produce no print to the screen.

An advantage of using semicolon is that the misspelling of restart will be obvious because it produces an output and a print to the screen, as in Restart;  or  retart; 

I wonder why you use `dsolve`instead of plain dsolve. There is no difference.

Doesn't the help page tell it as it is:

"The first integrals are set equal to generated global indexed variables K[i]  that denote arbitrary constants."

?

@tira I shall try again, but will begin with the conclusion:

Your expression V defined for y>0 has the limit 0 as y tends to zero from the right.

I do not use decimal numbers in order to show that the limit is indeed exactly zero.

restart:
t:=2/5; Pr:=71/100; Gamma:=1/5; Gr:=1/2;
## After this paste in your lines defining as you did all the terms v1, v2, ..., v8.
###
###Assuming that is done, continue with

C:=IntegrationTools:-Change(v4+v5-v3,u=y^2/x,[x]) assuming y>0; #Remember assuming
limit(C,y=0,right); #From the right in principle because the assumption y>0 was made
limit(v1-v2+C,y=0,right);
############Now we show that the 3 double integral terms v6, v7, v8 each has limit 0 as y->0:
## We shall use that
taylor(BesselI(1,x),x=0,4);
v6;
v6Int:=subs(int=Int,v6);
v6E:=simplify(eval(v6Int,BesselI=( ()->args[2]/2)));
# v6E is bounded above by
subs(y^2=0,v6E);
value(%);
v7;
v7Int:=subs(int=Int,-v7); #Notice minus: We are estimating absolute values
v7E:=simplify(eval(v7Int,BesselI=( ()->args[2]/2)));
# v7E is bounded above by
subs(y^2=0,v7E);
value(%);
v8;
v8Int:=subs(int=Int,v8);
v8E:=simplify(eval(v8Int,BesselI=( ()->args[2]/2)));
# v8E is bounded above by
subs(y^2=0,v8E);
value(%);
#Each of the upper bounds found for the absolute values of v6, v7, v8 has limit 0 as y->0.
##This finishes the proof of the conclusion stated in the beginning.

@tira The integrals in A-B of the same kind as v3 amount to v4+v5-v3.
Do the same substitution to that part of A-B :

IntegrationTools:-Change(v4+v5-v3,u=y^2/x,[x]) assuming y>0; ##Notice assumption
eval(%,y=0);
#Result 0.446573966
#  Then if you look at the contribution to A-B consisting of v1-v2, then you get
evalf(eval(v1-v2,y=0));
# except for roundoff likely the same.
#However, if y<0 is allowed in your setup, then you have a problem, because while the former result changes sign with y the latter does not.

Now left is the contribution to A-B from the double integrals v7+v8-v6. But because BesselI(1,x)  has a taylor expansion about 0 starting with x/2 each of the integrals itself is bounded in y thus with the factor y in front each of v6,v7,v8 tend to zero as y->0. 

Thus (except for the rounding error mentioned, which could be looked into) we see that the limit of A-B from the right as y->0 is 0.

The limit from the left is not zero though. But you may be assuming that y>=0 ?

@Kitonum It should be pointed out that your procedure Euler assumes that

1.  the derivative has been isolated on the left
2.  the dependent variable is called y
3.  the independent variable is called x

You have an h there, h=0.1. Does that mean that you wanted to solve the problem numerically?