Question: how this code works?

When a and b are from 0 to 1 it works!! but from -1 to -0.05 not i wonder why?

a1:=1.9:alpha:=0.2:
ode1 := diff( y(x), x$2 ) + x^(-1)* diff( y(x), x )  = 2*a1*(1+x^(2-alpha)/((2-alpha)^2*GAMMA(2-alpha))):
    bc1 := y(0)=subs(r=-0.05,a1/4*(r^2-1)), y(-1)=0:
a:=-1:b:=-0.05:
N := 20:h := (b-a)/N:
X := k -> a+k*h: 
'X'[k] = X(k):
Yp  := k -> (y[k+1]-y[k-1])/2/((1+alpha)*GAMMA(1+alpha)*h^(alpha)):Ypp := k -> (y[k+1]-2*y[k]+y[k-1])/((2-alpha)*GAMMA(2-alpha)*(h^(2-alpha))^2):
for k from 1 to N-1 do
eq[k] := eval( ode1,
                    {x=X(k), y(x)=y[k],
                     diff(y(x),x)=Yp(k),
                     diff(y(x),x$2)=Ypp(k)} ):
    end do:

eq[0] := y[0] = rhs(bc1[1]):
eq[N] := y[N] = rhs(bc1[2]):
    fd_sol1 := fsolve( {seq( eq[k], k=0..N )}, {seq( y[k], k=0..N )} ):
fd_table1 := eval( seq([X(k),y[k]],k=0..N), fd_sol1 ):Matrix([fd_table1]):
infolevel[dsolve] := 3:
#exact_sol1 := combine(dsolve( { ode1, bc1 }, y(x) )):
infolevel[dsolve] := 0:

P1:=plot([fd_table1], x=a..b,color=[black], linestyle = solid, symbolsize=16,axes=boxed);