In an expansion around the point \(z=1\), the recurrence relation is: \[a_{n+1}=\frac{1}{n+1}\, a_n\]
Sensemaking: The solution will look like \[y(z)= A\left[a_0+a_1 (z-1)+a_2 (z-1)^2+a_3 (z-1)^3+a_4 (z-1)^4+\dots\right]\] where my job is to find and plug in the \(a_n\)'s. The given recurrence relation involves only the indices \(n\) and \(n+1\) as would be typical for a first order ODE, so I expect one unknown constant, call it \(A\).
Using a recurrence relation, you will never find a value for the first one or two coefficents. These early coefficients will become the unknown constants in the solution of the homogeous differential equation. So let's set \(a_0=A\).
Then, the recurrence relation for \(n=0\) tells me that \[a_1=\frac{1}{0+1}\, a_0=a_0=A.\]
The recurrence relation for \(n=1\) tells me that \[a_2=\frac{1}{1+1}\, a_1= \frac{1}{2}\, a_1.\] But I can't stop there. I need every coefficient in terms of \(a_0=A\), so I plug in the result from \(n=0\): \[a_2=\frac{1}{2}\, a_1=\frac{1}{2}\, a_0=\frac{1}{2}\, A.\]
Continuing in this way, I find that:
\begin{align*} a_3 &= \frac{1}{3}\, a_2 =\frac{1}{3}\frac{1}{2}\, A =\frac{1}{3!}\, A\\ a_4 &= \frac{1}{4}\, a_3 =\frac{1}{4}\frac{1}{3!}\, A =\frac{1}{4!}\, A \end{align*}
Notice that I am leaving the coefficients with factorials in them so that I can see the pattern that is developing.
I now have 5 nonzero coefficients, so I can stop calculating, but I'm not done. I still need to write out all the terms in the power series. I'll call the unknown function that I am trying to find \(f\). I'm told that the recurrence relation is for an expansion around \(z=1\), so I expect to have powers of \(z-1\). Altogether: \[y(z)= A\left[1+(z-1)+\frac{1}{2!}\, (z-1)^2+\frac{1}{3!}\, (z-1)^3+\frac{1}{4!}\, (z-1)^4+\dots\right]\] Don't forget the \(a_0\) term which looks like \((z-1)^0=1\)! Note that in this case, because there is an \(n=0\) term, five nonzero terms means that this solution is correct to fourth order.
In an expansion around the point \(z=0\), the recurrence relation is: \[a_{n+2}=-\frac{(5-n)(6+n)}{(n+2)(n+1)}\, a_n\]
Sensemaking: The solution will look like \[y(z)= A\left[a_0+a_2 z^2+a_4 z^4+a_6 z^6+\dots\right] +B\left[a_1 z+a_3 z^3+a_5 z^5 +a_7 z^7+\dots\right]\] where my job is to find and plug in the \(a_n\)'s. The given recurrence relation involves only the indices \(n\) and \(n+2\) as would be typical for a second order ODE, so I expect two unknown constants, call them \(A\) and \(B\). One of the solutions will be even and one of the solutions will be odd.
Letting \(a_0=A\), where \(A\) is an unknown constant and using the same strategy as above, for \(n=0\), \(2\) and \(4\), I get:
\begin{align*} a_2 &= -\frac{5\cdot 6}{2\cdot 1}\, a_0= -\frac{5\cdot 6}{2\cdot 1}\, A \\ a_4 &= -\frac{3\cdot 8}{4\cdot 3}\, a_2 =\frac{5\cdot 3\cdot 6\cdot 8}{2\cdot 4\cdot 1\cdot 3}\, A \\ a_6 &= -\frac{1\cdot 10}{6\cdot 5}\, a_4 =-\frac{5\cdot 3\cdot 1\cdot 6\cdot 8\cdot 10}{2\cdot 4\cdot 6\cdot 1\cdot 3\cdot 5}\, A \\ \end{align*}
I will never get the odd terms from the sequence above. Instead, letting \(a_1=B\), where \(B\) is an unknown constant and using the same strategy as above, for \(n=1\), \(3\) and \(5\), I get:
\begin{align*} a_3 &= -\frac{4\cdot 7}{3\cdot 2}\, a_1= -\frac{4\cdot 7}{3\cdot 2}\, B \\ a_5 &= -\frac{2\cdot 9}{5\cdot 4}\, a_3 =\frac{4\cdot 2\cdot 7\cdot 9}{3\cdot 5\cdot 2\cdot 4}\, B \\ a_7 &= -\frac{0\cdot 11}{7\cdot 6}\, a_5 =0\, B \\ \end{align*}
Notice that the \(B\) series terminates. All the coefficients beyond \(a_5\) will be zero.
Putting it all together, I get: \begin{align*} y(z)&= A\left[1-\frac{5\cdot 6}{2\cdot 1} z^2+\frac{5\cdot 3\cdot 6\cdot 8}{2\cdot 4\cdot 1\cdot 3} z^4 -\frac{5\cdot 3\cdot 1\cdot 6\cdot 8\cdot 10}{2\cdot 4\cdot 6\cdot 1\cdot 3\cdot 5} z^6+\dots\right]\\ &+B\left[1 z-\frac{4\cdot 7}{3\cdot 2} z^3+\frac{4\cdot 2\cdot 7\cdot 9}{3\cdot 5\cdot 2\cdot 4} z^5 \right] \end{align*}
Notice that if you are asked for a certain number of nonzero terms and the series terminates, you should just give as many nonzero terms as there are.
Consider the differential equation \(y^{\prime\prime} - 2y^{\prime} + y = 0\).
Our answer at the end should have the form \(y(x) \approx A(c_0 + c_1x + c_2x^2 + c_3x^3 + c_4x^4 + c_5x^5) + B(d_0 + d_1x + d_2x^2 + d_3x^3 + d_4x^4 + d_5x^5)\). There are two solutions because it is a second-order differential equation. Our task is to find the \(c\) and \(d\) coefficients.
We start by assuming\begin{equation*} y(x) = \sum_{n=0}^\infty c_nx^n, \end{equation*}
giving
\begin{equation*} y^\prime(x) = \sum_{n=1}^\infty nc_nx^{n-1} \end{equation*}
and
\begin{equation*} y^{\prime\prime}(x) = \sum_{n=2}^\infty n(n-1)c_nx^{n-2} \end{equation*}
Rearranging the terms in the sums gives us:
\begin{equation*} \sum_{n=2}^\infty n(n-1)c_nx^{n-2} - 2\sum_{n=1}^\infty nc_nx^{n-1} + \sum_{n=0}^\infty c_nx^n = 0 \end{equation*} \begin{equation*} \sum_{n=0}^\infty ((n+2)(n+1)c_{n+2} - 2(n+1)c_{n+1} + c_n)(x^n) = 0 \end{equation*}
The parentheses must be equal to 0 for all \(n\), so:
\begin{equation*} c_{n+2} = \frac{2(n+1)c_{n+1} - c_n}{(n+1)(n+2)} \end{equation*}
We have the freedom to choose both \(c_0\) and \(c_1\). One option is to leave \(c_0\) and \(c_1\) as arbitrary constants, calculate the rest of the coefficients, then group like terms in order to get our two solutions. However, we can see from the recursion relation that \(c_{n+2}\) depends on both \(c_{n}\) and \(c_{n+1}\), so each coefficient will include a mix of \(c_0\) and \(c_1\). It's doable, but requires more algebra.
Another option is to consider \(c_0\) and \(c_1\) separately. That is, one solution can come from letting \(c_1=0\) and leaving \(c_0\) arbitrary and the second solution from letting \(c_0=0\) and \(c_1\) arbitary. Let's start by letting \(c_1=0\). This gives:
\begin{align*} &n=0: c_2 = -c_0/2 \\ &n=1: c_3 = 2c_2/3 = -c_0/3 \\ &n=2: c_4 = c_3/2 - c_2/12 = -c_0/8 \\ &n=3: c_5 = 2c_4/5 - c_3/20 = -c_0/30 \\ &n=4: c_6 = c_5/3 - c_4/30 = -c_0/144 \end{align*}
\(y_0(x) \approx c_0(1 - x^2/2 - x^3/3 - x^4/8 - x^5/30)\) does not immediately look like a familiar power series but it is still a valid solution. To get the second solution, we let \(c_0 = 0\) (and switch to \(d_n\) to remind ourselves that these coefficients are different from the first solution that we found):
\begin{align*} &n=0: d_2 = d_1 \\ &n=1: d_3 = d_1/2 \\ &n=2: d_4 = d_3/2 - d_2/12 = d_1/6 \\ &n=3: d_5 = 2d_4/5 - d_3/20 = d_1/24 \\ &n=4: d_6 = d_5/3 - d_4/30 = d_1/120 \\ \end{align*}
This looks very familiar! \(y_1(x) \approx d_1(x + x^2 + x^3/2 + x^4/6 + x^5/24) = d_1xe^x\).
The final answer puts these two solutions together as: \[y(x) \approx c_0(1 - x^2/2 - x^3/3 - x^4/8 - x^5/30) + d_1(x + x^2 + x^3/2 + x^4/6 + x^5/24)\]
As another option for a solution, we can try \(c_0=c_1\). Making this choice for the first solution might help you get a more recognizeable power series.
This equation can also be solved using the method of constant coefficients. The characteristic polynomial is: \(r^2 - 2r + 1 = (r - 1)^2\). Since this equation has one repeated solution, an additional independent solution must be added in the form of \(y(x) = Ae^x + Bxe^x\). The second term matches one of our solutions from before. The other term can be found from the series solution either by an appropriate superposition of the two solutions we found, or by going back to the series coefficients and choosing the initial condition of \(c_1 = c_0\).
In the previous problem, we found the normalized energy eigenstates \begin{align} \left|{n_x n_y}\right\rangle &\doteq \sqrt{\frac{2}{L_x}}\sqrt{\frac{2}{L_y}}\sin\left( \frac{n_x\pi x}{L_x} \right)\sin\left( \frac{n_yy\pi y}{L_y} \right)e^{-iE_{n_xn_y}t/\hbar}\\ E_{n_xn_y}&= \frac{\pi^2\hbar^2}{2m^2}\left( \frac{n_x^2}{L_x^2}+\frac{n_y^2}{L_y^2} \right) \end{align}
At this point we are finally prepared to write the full time-dependent wave function for the initial state \(\psi(x,y,0)\). Recall that we can do this by expanding the function in the energy basis, that is \begin{equation} \psi(x,y,t) = \sum_{n_x}\sum_{n_y}c_{n_xn_y}\sqrt{\frac{2}{L_x}}\sin\left( \frac{n_x\pi x}{L_x} \right) \sqrt{\frac{2}{L_y}} \sin\left( \frac{n_yy\pi y}{L_y} \right)e^{-iE_{n_xn_y}t/\hbar} \end{equation} which we would write in “ket-language” as \begin{equation} |\psi\rangle = \sum_{n_x}\sum_{n_y}c_{n_xn_y}e^{-iE_{n_xn_y}t/\hbar}|n_x,n_y\rangle \end{equation} All that we need to determine are the coefficients \(c_{n_xn_y}\). By now, this process should be familiar to you. The coefficients are given by \begin{align} c_{n_xn_y} &= \langle n_x,n_y|\psi \rangle\\ &= \int_0^{L_y}\int_0^{L_x}\sqrt{\frac{2}{L_x}}\sin\left( \frac{n_x\pi x}{L_x} \right) \sqrt{\frac{2}{L_y}} \sin\left( \frac{n_yy\pi y}{L_y} \right)\nonumber \\ &\qquad\times\frac{30}{\sqrt{L_xL_y}}(L_xx-x^2)(L_yy-y^2)\;dxdy \\ &= \frac{240\Big(-1+(-1)^{n_x}\Big)\Big(-1+(-1)^{n_y}\Big)}{n_x^3n_y^3\pi^6} \\ &= \begin{cases} \displaystyle{\frac{960}{n_x^3n_y^3\pi^6}}, & n_x,n_y\text{ are both odd} \\\ \;\;0, & \text{ otherwise} \end{cases} \end{align} It is perfectly reasonable to use Mathematica for this type of integral. If you do, take advantage of the Assumptions feature to tell Mathematica that \(n_x,n_y\) are integers (see attached code)
Our final, time-dependent wave function is given by \begin{equation} \psi(x,y,t) = \mathop{\sum\sum}_{n_x,n_y\text{odd}} \frac{960}{n_x^3n_y^3\pi^6}\sqrt{\frac{2}{L_x}}\sin\left( \frac{n_x\pi x}{L_x} \right) \sqrt{\frac{2}{L_y}} \sin\left( \frac{n_yy\pi y}{L_y} \right)e^{-iE_{n_xn_y}t/\hbar} \end{equation}