Series Solutions

Series solutions near an ordinary point for a second order, linear, homogenous DE
    (1-x) y'' + y = 0

Clear[y,x]
DSolve[(1-x) y''[x]+y[x]==0,y[x],x]

what can we tell from the DE about a solution?

Note that in standard form, the DE is y'' + (0) y' + ((1/(1-x)))y = 0.  The coefficients are p(x) = 0, q(x) = (1/(1-x)), g(x) =0.  These are analytic (i.e.,Taylor Series with a positive radius of convergence, i.e., have an infinite number of continuous derivatives) at every point of the real line except (Subscript[x, 0])=1.   In other words,  (Subscript[x, 0])=1 is a singular point of the DE and every other point of the real line is an ordinary point.  In particular, (Subscript[x, 0])=0 is an ordinary point, and we may seek a solution of the form y(x) = (∑_(n=0)^∞)(Subscript[a, n]) (x^n).  Theory only guarantees that such a series solution converges on the interval (-1,1).  The solution may be good outside of this interval.

3 terms in each of 2 linearly independent solns

First state the form of the solution to a desired number of terms.

Clear[y,a,y1,y2,i,x,lhs, ysoln]

y=Sum[a[i] x^i, {i,0,6}] + O[x]^7

Calculate the necessary derivatives

y1=D[y,x]

y2=D[y1,x]

Substitute y and its derivatives into the differential equation.

lhs=Expand[(1-x)*y2+y]

Equate the coefficients for corresponding powers of x.

LogicalExpand[lhs==0]

Since the DE is second order linear homogeneous, we should be able to solve for all but the first two terms in the series solution.

Solve[%,{a[2],a[3],a[4],a[5],a[6]}]

Substitute these back into y.

{ysoln}=y/.%

Collect the terms for the two linearly independent series solutions.

ysoln=Collect[ysoln,{a[0],a[1]}]

ysoln

Look at the portions of this solution.  The linearly independent solutions to the homogenous equation are 
y1(x) = 1-((x^2)/2)-((x^3)/6)-((x^4)/24)-((x^5)/60)-((7 (x^6))/720)+. . .  and y2(x) = x-((x^3)/6)-((x^4)/12)-((x^5)/24)-((x^6)/40)+ . . .

some of these partial sums

ysoln/.{x->0}

What is y' (0)?

D[ysoln,x]

In this case, a[0] is y(0), the y-intercept, and a[1] is the slope of the curve at that point.

D[ysoln,x]/.{x->0}

Do not worry about the code here; we generate some initial values and some solution curves.  Here we have chosen some points with different y-intercepts and nonnegative slopes at the initial value of x=0.

Clear[pts,ptList,fcnList,c1,c2];
ptList={};fcnList={};

Do[AppendTo[ptList,{0,c1}],{c1,-4,4,2}]
For[c1=-4,c1≤4,
For[c2=0,c2≤3,c2++,AppendTo[fcnList,(1-((x^2)/2)-((x^3)/6)-((x^4)/24)-((x^5)/60)-((7 (x^6))/720)) c1+(x-((x^3)/6)-((x^4)/12)-((x^5)/24)-((x^6)/40)) c2]];
c1=c1+2]
ptList
fcnList

pts =  Graphics[{PointSize[Large],Point[ptList]}];
Show[{Plot[fcnList,{x,-3,4},PlotRange->{-9,11}],pts}]

Aiding in handwork

y(x) = (∑_(n=0)^∞)(Subscript[a, n]) (x^n).  y'(x)  = (∑_(n=1)^∞)n (Subscript[a, n]) (x^(n-1)).   y''(x)  = (∑_(n=2)^∞)n(n-1)(Subscript[a, n]) (x^(n-2)). 
0 = (1 - x) y''(x) + y(x)  = (∑_(n=2)^∞)n(n-1)(Subscript[a, n]) (x^(n-2)) - x (∑_(n=2)^∞)n(n-1)(Subscript[a, n]) (x^(n-2)) + (∑_(n=0)^∞)(Subscript[a, n]) (x^n) 
=  (∑_(n=2)^∞)n(n-1)(Subscript[a, n]) (x^(n-2)) - (∑_(n=2)^∞)n(n-1)(Subscript[a, n]) (x^(n-1)) + (∑_(n=0)^∞)(Subscript[a, n]) (x^n)  (let j=n-2 in series 1, j=n-1 in series 2, j=n in series 3)
= (∑_(j=0)^∞)(j+2)(j+1)(Subscript[a, j+2]) (x^j) - (∑_(j=-0)^∞)(j+1)j (Subscript[a, j+1]) (x^j) + (∑_(j=0)^∞)(Subscript[a, j]) (x^j)
= (∑_(j=0 )^∞){  (j+2)(j+1)(Subscript[a, j+2]) - (j+1)j (Subscript[a, j+1]) + (Subscript[a, j])}(x^j)

Table[a[j+2] ==j/(j+2) * a[j+1] - 1/((j+1)*(j+2)) a[k],{j,1,4}]

a[3]==-((a[1])/6)+((a[2])/3)/.a[2]->-1/2 a[0]

a[4]==-((a[2])/12)+((a[3])/2) /.{a[3]->- ((a[0])/6)-((a[1])/6),a[2]-> -a[0]/2}

Simplify[%]

Solve[%,a[4]]

a[5]==-((a[3])/20)+((3 a[4])/5)/.{a[3]-> -((a[0])/6)-((a[1])/6),a[4]-> (1/24) (-a[0]-2 a[1])}

Simplify[%]

Solve[%,a[5]]

a[6]==-((a[4])/30)+((2 a[5])/3)/.{a[5]->(1/120) (-2 a[0]-5 a[1]),a[4]-> (1/24) (-a[0]-2 a[1])}

Simplify[%]

Solve[%,a[6]]

10 terms in each solution

Of course, this same technique is used if we wish to take these to 10 terms each.

Clear[y,a,y1,y2,i,x,lhs, ysoln, eqnset, varset, u1, u2]

y=Sum[a[i] x^i, {i,0,20}] + O[x]^21

y1=D[y,x]

y2=D[y1,x]

lhs=Expand[(1-x)*y2+y]

eqnset=LogicalExpand[lhs==0]

varset=Table[a[i],{i,2,20}]

Solve[eqnset,varset]

{ysoln}=y/.%

ysoln=Collect[ysoln,{a[0],a[1]}]

u1=Coefficient[ysoln,a[0]]

u2=Coefficient[ysoln,a[1]]

The two linearly independent solutions to our homogeneous DE are u1(x) + ... and u2(x) + ...

some of these partial sums

We generate some initial values and some solution curves.

Clear[pts,ptList,fcnList,c1,c2];
ptList={};fcnList={};

Do[AppendTo[ptList,{0,c1}],{c1,-4,4,2}]
For[c1=-4,c1≤4,
For[c2=0,c2≤3,c2++,AppendTo[fcnList,u1 c1+u2 c2]];
c1=c1+2]

pts =  Graphics[{PointSize[Large],Point[ptList]}];
Show[{Plot[fcnList,{x,-3,4},PlotRange->{-9,11}],pts}]

Series solutions near an ordinary point for a second order, linear, nonhomogenous DE
    (1-x) y'' + y = 3 (x^2) - x + 4

what can we tell from the DE about a solution?

Note that in standard form, the DE is y'' + (0) y' + ((1/(1-x)))y = (((3(x^2)-x+4)/(1-x))).  The coefficients are p(x) = 0, q(x) = (1/(1-x)), g(x) =((3(x^2)-x+4)/(1-x)).  These are analytic (i.e.,Taylor Series with a positive radius of convergence, i.e., have an infinite number of continuous derivatives) at every point of the real line except (Subscript[x, 0])=1.   In other words,  (Subscript[x, 0])=1 is a singular point of the DE and every other point of the real line is an ordinary point.  In particular, (Subscript[x, 0])=0 is an ordinary point, and we may seek a solution of the form y(x) = (∑_(n=0)^∞)(Subscript[a, n]) (x^n).  Theory only guarantees that such a series solution converges on the interval (-1,1).  The solution may be good outside of this interval.

3 terms in each of 2 linearly independent solns

First state the form of the solution to a desired number of terms.

Clear[y,a,y1,y2,i,x,lhs, ysoln]

y=Sum[a[i] x^i, {i,0,6}] + O[x]^7

Calculate the necessary derivatives

y1=D[y,x]

y2=D[y1,x]

Substitute y and its derivatives into the differential equation.

lhs=Expand[(1-x)*y2+y]

Equate the coefficients for corresponding powers of x.

LogicalExpand[lhs == 3 x^2 - x + 4]

Since the DE is second order, we should be able to solve for all but two terms in the series solution.  Since the DE is nonhomogeneous, we are not guaranteed to find all in terms of the first two, but it works here.

Solve[%,{a[2],a[3],a[4],a[5],a[6]}]

Substitute these back into y.

{ysoln}=y/.%

Collect the terms for the two linearly independent series solutions; we're only showing a certain number of terms in each.

ysoln=Collect[ysoln,{1, a[0],a[1]}]

Look at the portions of this solution.  The linearly independent solutions to the associated homogenous equation are 
y1(x) = 1-((x^2)/2)-((x^3)/6)-((x^4)/24)-((x^5)/60)-((7 (x^6))/720)+. . .  and y2(x) = x-((x^3)/6)-((x^4)/12)-((x^5)/24)-((x^6)/40)+ . . . ;

the particular solution is yp(x) = 2 (x^2)+((x^3)/2)+((x^4)/3)+((7 (x^5))/40)+((19 (x^6))/180)+...

u1=Coefficient[ysoln,a[0]]

u2=Coefficient[ysoln,a[1]]

yp = ysoln - u1 a[0] - u2 a[1]

some of these partial sums

ysoln/.{x->0}

What is y' (0)?

D[ysoln,x]

In this case, a[0] is y(0), the y-intercept, and a[1] is the slope of the curve at that point.

D[ysoln,x]/.{x->0}

Do not worry about the code here; we generate some initial values and some solution curves.  Here we have chosen some points with different y-intercepts and nonnegative slopes at the initial value of x=0.

Clear[pts,ptList,fcnList,c1,c2];
ptList={};fcnList={};

Do[AppendTo[ptList,{0,c1}],{c1,-4,4,2}]
For[c1=-4,c1≤4,
For[c2=0,c2≤3,c2++,AppendTo[fcnList,u1 c1+u2 c2 + yp] ];
c1=c1+2]
ptList
fcnList

pts =  Graphics[{PointSize[Large],Point[ptList]}];
Show[{Plot[fcnList,{x,-3,4},PlotRange->{-9,11}],pts}]

Aiding in handwork

y(x) = (∑_(n=0)^∞)(Subscript[a, n]) (x^n).  y'(x)  = (∑_(n=1)^∞)n (Subscript[a, n]) (x^(n-1)).   y''(x)  = (∑_(n=2)^∞)n(n-1)(Subscript[a, n]) (x^(n-2)). 
3 (x^2) - x + 4 = (1 - x) y''(x) + y(x)  = (∑_(n=2)^∞)n(n-1)(Subscript[a, n]) (x^(n-2)) - x (∑_(n=2)^∞)n(n-1)(Subscript[a, n]) (x^(n-2)) + (∑_(n=0)^∞)(Subscript[a, n]) (x^n) 
=  (∑_(n=2)^∞)n(n-1)(Subscript[a, n]) (x^(n-2)) - (∑_(n=2)^∞)n(n-1)(Subscript[a, n]) (x^(n-1)) + (∑_(n=0)^∞)(Subscript[a, n]) (x^n)  (let j=n-2 in series 1, j=n-1 in series 2, j=n in series 3)
= (∑_(j=0)^∞)(j+2)(j+1)(Subscript[a, j+2]) (x^j) - (∑_(j=-0)^∞)(j+1)j (Subscript[a, j+1]) (x^j) + (∑_(j=0)^∞)(Subscript[a, j]) (x^j)
= (∑_(j=0 )^∞){  (j+2)(j+1)(Subscript[a, j+2]) - (j+1)j (Subscript[a, j+1]) + (Subscript[a, j])}(x^j)

coefficient of x^0

Solve[4 == 2 a[2] + a[0],a[2]]

Expand[%]

subs1 = {a[2]->2-((a[0])/2)}

(a[2+j](1+j) (2+j) -j (1+j) a[1+j]+a[j])/. {j->1}

coefficient of x^1

Solve[-1 == a[1]-2 a[2]+6 a[3],a[3]]

Expand[%]

%/.subs1

Expand[%]

subs2 = {a[2]->2-((a[0])/2),a[3]->(1/2)-((a[0])/6)-((a[1])/6)}

(a[2+j](1+j) (2+j) -j (1+j) a[1+j]+a[j])/. {j->2}

coefficient of x^2

Solve[3 == a[2]-6 a[3]+12 a[4],a[4]]

%/.subs2

Expand[%]

subs3 = {a[2]->2-((a[0])/2),a[3]->(1/2)-((a[0])/6)-((a[1])/6),a[4]->(1/3)-((a[0])/24)-((a[1])/12)}

coefficient of x^n  for n>=3

Solve[a[2+j](1+j) (2+j) -j (1+j) a[1+j]+a[j] == 0, a[2+j]]

Simplify[%]

Factor[%]

{a[2+j]->((-a[j]+j(j+1) a[1+j])/((1+j) (2+j)))}/. {j->3}

Expand[%]

% /.subs3

Expand[%]

subs4 = {a[2]->2-((a[0])/2),a[3]->(1/2)-((a[0])/6)-((a[1])/6),a[4]->(1/3)-((a[0])/24)-((a[1])/12),a[5]->(7/40)-((a[0])/60)-((a[1])/24)}

{a[2+j]->((-a[j]+j(j+1) a[1+j])/((1+j) (2+j)))}/. {j->4}

Expand[%]

% /.subs4

Expand[%]