Differential Equation Solution By Power Series 
Solve $(1 + x)y' = py;\ \ \ y(0) = 1$, where $p$ is an arbitrary constant.

First I plugged in the guess $y = \sum_{n = 0}^\infty a_n x^n$:
$(1 + x)(\sum_{n = 0}^\infty a_n x^n)' = p\sum_{n = 0}^\infty a_n x^n$
Then I expanded the derivative and multiplication:
$\sum_{n = 0}^\infty n a_n x^{n - 1} + \sum_{n = 0}^\infty n a_n x^n = p\sum_{n = 0}^\infty a_n x^n$
Then I shifted the left index (the first term yielding $0$ allows the lower bound to remain $0$) and algebraically combined the summations:
$\sum_{n = 0}^\infty (n + 1)a_{n + 1} x^n + (n - p)a_n x^n = 0$
This leads to the following recurrence relation:
$a_{n + 1} = \frac{p - n}{n + 1}a_n$
Thus for various values of $n$:
$a_1 = p a_0$, $a_2 = \frac{p(p - 1)}{2}a_0$, $a_3 = \frac{p(p - 1)(p - 2)}{6} a_0$, etc.
So applying definitions for the exponential taylor series and falling factorial, the guessed solution would be:
$y = \sum_{n = 0}^\infty \frac{p! a_0 x^n}{n! (p - n)!} = \sum_{n = 0}^\infty a_0 e^x p^{\underline n}$
Solving the initial value problem:
$1 = \sum_{n = 0}^\infty a_0 e^0 p^{\underline n} \implies a_0 = \frac{1}{\sum_{n = 0}^\infty p^{\underline n}}$
My final solution is:
$y = \frac{\sum_{n = 0}^\infty e^x p^{\underline n}}{\sum_{n = 0}^\infty p^{\underline n}}$
However, the answer is supposed to be $y = (1 + x)^p$.  Are these identical, or did I make an error somewhere?
 A: Your solution
$$
y = \sum_{n = 0}^\infty \frac{p! a_0 x^n}{n! (p - n)!} 
$$
is mostly the binomial series for $a_0\,(1+x)^p$, assuming you know what you are doing (it looks like you are treatng $p$ as an integer, which it isn't necessarily according to the question as you posted it).
The "equality"
$$
\sum_{n = 0}^\infty \frac{p! a_0 x^n}{n! (p - n)!}= \sum_{n = 0}^\infty a_0 e^x p^{\underline  n}
$$
is not true; it's missing the $n!$ in the denominator on the right.
A: Hint: You know that the solution (by separation of variables) is:
$y(x)=(x+1)^p$
Now try to use the binomial theorem, expand the above, and equate coefficients. You will see that:
$a_0=1$
$a_1=pa_0=p= \binom{p}{1}$
$a_2=\frac{p(p-1)}{2}a_0= \binom{p}{2}$
.
.
.
$a_k=\binom{p}{k}$
As the others have pointed out, $a_k=0$ for $k\geq p+1$. This is because the solution is a polynomial.
So your recursion is correct, but after you applied the exponential Taylor series, things took a different turn.
Moral of the story: If you have an easier method, use it to double check your work :)
Edit: The above assumes that $p$ is an integer. If $p$ is an arbitrary scalar value (i.e. real, or even complex) it turns out we can in fact use a similar formula, we just have to be careful of how we are computing the binomial coefficients.
Indeed, we have that for arbitrary $p$,
$$y(x)= \sum_{k=0}^{\infty}\binom{p}{k}x^k$$ for $|x|<1$ and $x,p \in \mathbb{R}$
I did a little investigating, and found this to be helpful:
Generalized Binomial Expansion of $\left(1+x \right )^{y}$
Edit: Looks like there are multiple formulas for the "generalized binomial coefficients", but digging into Wolfram Mathworld, I came across:
$$\binom{p}{k}:= \frac{(p)_k}{k!}=\frac{p(p-1)(p-2)\cdots(p-k+1)}{k(k-1)\cdots 2 \cdot 1}$$
Look at formulas $(3),(4)$ in:
https://mathworld.wolfram.com/BinomialCoefficient.html
