I solve such questions using the Carleman-matrix-concept (which you also find mentioned in the wikipedia article). Carleman-matrices contain in its rows the coefficients of a function in its power series representation, and of its powers. So the Carleman-matrix, say C contains the coefficients of $f(x)^0$ (which is the constant 1), $f(x)$, $f(x)^2$ ... in its rows.
In your case where $f(x)=1 - 1 x$ we have $$ C= \begin{bmatrix}
1&.&.&. \\
1&-1&.&. \\
1&-2&1&. \\
1&-3&3&1 \\
\end{bmatrix} \qquad \text{ for } \qquad \begin{array} {}
f(x)^0 &=&1 \\
f(x)&=&1-x \\
f(x)^2&=&1-2x+x^2 \\
f(x)^3&=&1-3x+3x^2-x^3 \\
\end{array} $$
such that with a "vandermonde"-type column-vector $V(x) = [1,x,x^2,x^3,...]$ of appropriate size we shall have
$$ C \cdot V(x) = V(f(x)) $$
Then if you find a diagonalization of C such that $M^{-1}\cdot D \cdot M = C $ where $D$ is diagonal and $M$ and $M^{-1}$ are triangular, then $M$ in its second row contains the coefficient of the Schröder-function and that in $M^{-1}$ in its second row its inverse. In our case we find that
$$ M^{-1} =\begin{bmatrix}
1 & . & . & . \\
1/2 & 1/2 & . & . \\
1/6 & 1/2 & 1/3 & . \\
0 & 1/4 & 1/2 & 1/4
\end{bmatrix} $$
, with $D=\operatorname{diag}([1,-1,1,-1])$ and
$$ M = \begin{bmatrix}
1 & . & . & . \\
-1 & 2 & . & . \\
1 & -3 & 3 & . \\
-1 & 4 & -6 & 4
\end{bmatrix}$$
is a possible solution.
(Here $M^{-1}$ can be recognized as the set of coefficients of integrals of the bernoulli-polynomials when we extend the dimension of the matrix infinitely)
In general, the eigenvector-matrix $M$ can be understood as limit of the n'th power of $C$ scaled by the reciprocal of the n'th power of $f'(0)$ when n goes to infinity, and so the Schröder-function as limit of the n'th iterate of $f(x)$ divided by $f'(0)^n$ where $n \to \infty \qquad $ - but can furtherly be scaled by an arbitrary constant factor $\gamma \ne 0$
Remark: Your example which requires only matrix size of $n \times n= 2 \times 2$ is much easier, but then one wouldn't see the general principle (and the relation to the bernoulli-polynomials, so I used the bigger matrix size with n=4 here.