I'm confused as to why dividing $p(x) = q(x) / q(x_1)$ works. I understand that $q(x)$ is designed such that any $x_i$ for $2 \leq i \leq d+1$ will make $q(x)=0$. I know that $q(x_1)$ will not make $q(x) =0$ because there is no $(x-x_1)$. So, you keep generating a $q(x)$ that is non-zero for $x_i$ for all $2 \leq i \leq d+1$, and then divide by their respective $q(x_i)$ ? How does this give us $p(x)$ ?
The method described in your question explains how to obtain a solution when $y_1=1$ and all other $y_i$ are zero, but of course it doesn't give you the general solution. But once you have this, there isn't much left to do.
Denote $p_1(x)$ the polynomial obtained by this method. Using the same construction, you can obtain other polynomials $p_i(x)$, which are such that $p_i(x_i)=1$ and $p_i(x_j)=0$ for all $i,j$.
Once you have all these polynomials, you can find the interpolating polynomial for any given $y_1,\dots y_{d+1}$; it is simply equal to $$p(x)=\sum_{k=1}^{d+1}y_kp_k(x).$$ Indeed, you can check that this polynomial is such that $$p(x_i)=\sum_{k=1}^{d+1}y_kp_k(x_i)=y_ip(x_i)=y_i$$ (the other terms in the sum being all zero).