Linear differential equation
Expansion at irregular singular point

Expansion around x  =  
      0  for   y''[x] - (1/(x^3)) y[x] = 0. Example 2 in Bender & Orszag Section 3.4.

Initial change of variable  y[x] = (e^s)[x]

If we replace y[x]  by  (e^(s[x]  : ))

y[x_]:=Exp[s[x]]

the general second order equation in y[x] (only LHS shown; RHS = 0)  gets transformed into one in s[x]

eq1=y''[x]+p[x] y'[x]+q[x] y[x]

\(TraditionalForm
    \`Multiplying\ through\ by\ e\^\(-s[x]\)\ cleans\ this\ \ up\ a\ 
      \(bit : \)\)

eq2=Expand[eq1*Exp[-s[x]]]

In our particular case,  p[x] and q[x] are defined as given below. The ODE for s[x] then simplifies:

p[x_]:=0 ; q[x_]:=-x^(-3) ; eq3=eq2

First approximation to s[x]

At this stage, the s''[x]-term is usually asymptotically small compated to the other two terms. We assume this to be the case, and solve the simplified first order ODE this leads to

DSolve[-(1/(x^3))+((s'[x])^2)==0,s[x],x]

The second order equation we started with has two possible solutions. We choose to consider the one that grows as x->0, i.e. choose

s[x_] := (2/(Sqrt[x]))

To check that our assumption about s''[x] small compared to S'[x]^2 and x^(-3), we compare these

{s''[x], s[x]^2, x^(-3)}

As x->0, the assumption was indeed OK.

Refinement of approximation to s[x]

Let c[x] be the first correction to s[x], i.e.

s[x_] := (2/(Sqrt[x])) + c[x]

and the equation for c[x] (still, we work with the LHS only;  RHS=0) becomes

eq3

which simplifies to

eq4=Expand[%]

\(TraditionalForm
    \`We\ now\ know\ that\ c[x] << \(x\^\(\(-1\)/2\)\)\ \(( ∼ s[x])\), 
    \ \ \ implying\ \ \(c'\)[x] << \(x\^\(\(-3\)/2\)\), \ \ 
    \(\(c'\)'\)[x] << \(x\^\(\(-5\)/2\)\).\ The\ last\ two\ terms\ in\ the\ 
      equation\ are\ small\ compared\ to\ the\ first\ one.\ The\ dominant\ 
      balance\ must\ occur\ between\ the\ first\ two\ terms.\ Hence\ we\ solve
      \)

DSolve[(3/(2 (x^(5/2))))-((2 c'[x])/((x^(3/2))))==0,c[x],x]

For next approximation to s[x], we therefore try

s[x_] := (2/(Sqrt[x])) + 3/4 Log[x] + d[x]

where d[x] will have to satisfy

Expand[eq3]

This time,  we know d[x] << Log[x],  
    implying d'[x] << (x^(-1)),  
    d''[x] << (x^(-2)). Again,  
    the dominant balance has to be between the first two terms,  so

DSolve[-(3/(16 (x^2)))-((2 d'[x])/((x^(3/2))))==0,d[x],x]

At this point, the corrections go to zero as x->0, i.e. we can no longer ignore the integration constant. We could continue in the same style, and get more corrections to the finction s[x], but it is now often more practical to aim instead for a refinemt of the original function y[x]. We have so far found that
y[x] ∼ Exp[(2/(Sqrt[x])) + 3/4 Log[x]+c[1]-((3 (Sqrt[x]))/16)] ∼ c[2]  (x^(3/4 ))(e^(2/(Sqrt[x]) ))(1 - ((3 (Sqrt[x]))/16)+...)
Next tast is to generate some more terms in the last expansion shown.

Refinement of approximation to y[x]

Substitution of

y[x_]  := 
      x^(3/4) Exp[2/x^(1/2)] (1 - 3/16*x^(1/2) + e[x])

into the the ODE for y[x] gives

Simplify[ y''[x] + p[x] y'[x] + q[x] y[x] ]

eq5 = Simplify[256 (x^(5/4)) Exp[-2/(Sqrt[x]) ]*%]

or simpler

With e[x] << (x^(1/2)),  
    e'[x] << (x^(-1/2)),  
    e''[x] << (x^(-3/2)),  
    the only term above that can reach the size of the first one 
      is the third 
      one  :    128 (-4 (Sqrt[x])) e'[x].  Hence,  
    we solve

DSolve[-15 (Sqrt[x])+128 (-4 (Sqrt[x])) e'[x]==0,e[x],x]

\(TraditionalForm
    \`\(The\ expansion\ for\ y[x]\ \((following\ the\ exponential\ factor)\)\ 
      is\ \(therefore(1\  - \ \(3\ \@x\)\/16 - \(15\ x\)\/512 + ...)\).\ It\ 
      looks\ very\ likely\ that\ it\ will\ continue\ as\ a\ series\ in\ 
      x\^\(1/2\).\ Calling\ this\ function\ w[x], \ 
    it\ will\ have\ to\ cause\ y[x]\ to\ satisfy\ \ \ \)\)

y[x_]  := x^(3/4) Exp[2/x^(1/2)] w[x]

i.e. the following expression needs to be zero

y''[x] + p[x] y'[x] + q[x] y[x]

This can be cleaned up a bit

eq7 = Expand[ % * Exp[-2/(Sqrt[x])]]

We can now try to get the first  n  (say, n=6) coefficients through

n = 6  ;   w[x_] = 1 + Sum[a[k] x^(k/2), {k, 1, n}]

Next, plug this into the equation for w[x]

Expand[eq7]

After getting rid of the leading factor (and throwing in an extra (x^(1/2)) so that we can put to zero in turn the coefficients for (x^(1/2)),  (x^(2/2)), (x^(3/2))  etc;  Maple's routine Ccoefficient does not like to extract the coefficient for (x^0))

eq6 = Expand[% * 16 (x^(7/4))]

we get the successive coefficients through

Do [ 
      a[i] = a[i] /. 
          Solve[Coefficient [eq6, x^(1/2), i] == 0, a[i]][[1]]; 
       Print[a[i]],  {i, 1, n}]

This expansion has apparently not encountered any problems. We can easily obtain any number of terms we wish. The algebra in this very last step (running out many coefficients) is trivial once we have the simple ODE that w[x] satisfies. In Bender & Orszag, pages 85-87, it is shown how we can do the recusrions so that we get all the coefficients in closed form.