How to make this polynomial the zero polynomial?(recursively)? Given a fixed $\beta \in \mathbb{R}$, I want to find the $c_0,...,c_n$ for arbitrary $n \in \mathbb{N}$ such that the polynomial 
\begin{align}P_n(z):=z(1-z) \left(n(n-1)z^{n-2}+\sum_{i=2}^{n-1}i(i-1)c_i z^{i-2} \right)+\left(\frac12+(4\beta-1)z-4\beta z^2\right)\left(nz^{n-1}+\sum_{i=1}^{n-1}ic_i z^{i-1} \right)+\\
+\left((2n+1) \beta(2z-1)+\beta^2(1-4z) -2 \beta(1-2\beta)z+\beta+c_n\right)\left(z^n+\sum_{i=0}^{n-1}c_i z^{i} \right)=0. 
\end{align}
Is there a way to find recursive or explicite representations of $c_0,...,c_n$?
The problem is apparently, that $c_n$ makes this thing non-linear, otherwise this would be just solving a linear system. 
 A: Condensed formulation
\begin{align}P_n(z):=z(1-z) p''(z)+\left(\tfrac12+(4β-1)z-4β z^2\right)p'(z)+
\left(4nβz+β^2-2nβ+c_n\right)p(z)=0. 
\end{align}
where
$p(z)=z^n+\sum_{i=0}^{n-1}c_iz^i$
so that the formally highest degree $n+1$ of this construction cancels in its coefficients. The overall structure of this equation is
$$
L_n(p)=-c_n\,p,
$$
a classical eigenvalue problem. Any solution $p(z)=a_nz^n+...+a_1z+a_0$ can be transformed into the required form by dividing out $a_n$, so $c_k=a_k/a_n$, $k=0,1,...,n-1$.

Example/Regression test:
$n=1$, $β=1$ as in the addendum to the question. Then, using simple arithmetics,
$$
L_1(z+c_0)
=0+\left(\tfrac12+3z-4 z^2\right)\cdot 1+(4z-1)(z+c_0)
=(2+4c_0)z+\tfrac12-c_0
$$
And thus $L_1p=-c_1p$ iff $2(1+2c_0)=-c_1$ and $\tfrac12-c_0=-c_1c_0=2c_0+4c_0^2$ or
$$
c_0=-\frac18 \left(3\pm\sqrt{17}\right)\in\{-0.89038820...,0.14038820320...\}
$$
In matrix form, one gets from $L_1(z)=2z+\tfrac12$ and $L_1(1)4z-1$ the system matrix
$$
\begin{bmatrix}
-1&0.5\\4&2
\end{bmatrix}
$$
with the eigenvalue equation $0=λ^2-λ-4=(λ-\tfrac12)^2-\tfrac{17}{4}$ resulting in
$$
c_1=-λ=-\tfrac12(1\pm\sqrt{17})
$$
and $4a_0+(2+c_1)a_1=0$ giving the same values for $c_0=a_0/a_1$ as above.

The eigenproblem in polynomial coefficients
Applied to the monomial basis, the linear differential operator gives
\begin{align}
L_n(z^k)&=k(k-1)(z^{k-1}-z^k)+k\left(\tfrac12z^{k-1}+(4β-1)z^k-4β z^{k+1}\right)+\left(4nβz^{k+1}+(β^2-2nβ)z^k\right)\\
&=(k^2-\tfrac k2)z^{k-1}+(-k^2-2β(n-2k)+β^2))z^k+4β(n-k)z^{k+1}
\end{align}
which results in a tridiagonal coefficient matrix.
For $n=1=β$, this gives the same coefficients as above.

Old idea, construct system from polynomial values
Assume that, contrary to the computations above, the coefficients of the linear operator are so complicated that it hast to be treated as essentially a black box. After setting the parameters $β,c_0,...,c_{n-1},c_n$, only evaluation is easily achievable.
Essentially, the system $(L_n-c_n)(p)=0$ is semi-linear in the sense that it is of the form $$A(c_n)\cdot c = b(c_n),$$
$c=(c_0,...,c_{n-1})$. The system matrix is of dimension $(n+1)\times n$ since the degrees of the actual output range from $0$ to $n$. For the system to be solvable the variable $c_n$ has to be chosen so that $rank(A(c_n),b(c_n))=n$ which gives a number $(n+1)$ of determinant equations. There are $n+1$ conditions for $n+1$ variables. Since $c_n$ occurs linearly in the coefficients, the system reduces to a generalized eigenvalue equation
$$
((1-c_n)[A(0),b(0)]+c_n[A(1),b(1)])\begin{bmatrix}
c_0\\...\\c_{n-1}\\1
\end{bmatrix}=0
$$
This is again solvable using standard algorithms. The matrices $A(0),b(0),A(1),b(1)$ are constructed as above by combining $c_n\in\{0,1\}$ with $(c_0,c_1,...,c_{n-1})\in\{e_0,...,e_{n-1}\}\subset\Bbb R^n$ and evaluating for a set of $n+1$ sampling points $z\in\{x_0,x_1,...,x_n\}$.

Introducing a bigger notation $P_n(c,c_n,z)$ for the polynomial, the matrix coefficients are
$$
A_{i,j}(c_n)=P_n(e_j,c_n,x_i)\text{ and }b_i(c_n)=P(0,c_n,x_i), i=0,1,...,n, \;j=0,1,...,n-1.
$$
which are linear in $c_n$.
The equation system for the polynomial then reads as
$$
\left(\sum_{j=0}^{n-1} A_{i,j}(0)c_j+b_i(0)\right)=-c_n\left(\sum_{j=0}^{n-1} (A_{i,j}(1)- A_{i,j}(0))c_j+(b_i(1)-b_i(0))\right)
$$
which is a generalized eigenvalue problem with $(c_0,...,c_{n-1},1)$ as eigenvector. If $[ΔA,Δb]=[A(1)-A(0),b(1)-b(0)]$ is invertible, it is equivalent to a normal eigenvalue problem, but even if that is not the case, the generalized problem can have solutions.
