The question asks to find the MLE of theta, for the following function:

This was my attempt:
And this was the solution:
Where did I go wrong in my working?
On
$L = (1 + \theta)^n \prod\limits_{i=1}^n x_{i}^{\theta} \Rightarrow \log L = n \log(1 + \theta) + \theta \sum\limits_{i=1}^n\log x_{i} \to \max\limits{\theta}$
$\log L^{\prime}_{\theta} = \frac{n}{1 + \theta} + \sum\limits_{i=1}^n \log x_{i} = 0 \Rightarrow \boxed{\hat{\theta} = \frac{-n}{\sum\limits_{i=1}^n \log x_{i}} - 1}$
The joint density is indexed by the observation number $i$, so you cannot write the step in red: $$\prod_{i=1}^n f(x_i \mid \theta) = \prod_{i=1}^n (1+\theta) x_i^\theta \color{red}{\ne (1+\theta)^n x^{n\theta}}.$$ This is because $x_1, x_2, \ldots, x_n$ are not necessarily the same number; they represent a realization of your sample of observed data from this distribution. Consequently, $$\prod_{i=1}^n (1+\theta) x_i^\theta = (1+\theta)^n \left(\prod_{i=1}^n x_i\right)^\theta.$$ The log-likelihood is therefore $$\ell(\theta \mid x_1, \ldots, x_n) = n \log (1 + \theta) + \theta \log \left(\prod_{i=1}^n x_i \right)= n \log (1 + \theta) + \theta \sum_{i=1}^n \log x_i.$$ In the last step, we used a generalized form of the property of logarithms $$\log ab = \log a + \log b, \quad a, b > 0.$$ The rest is straightforward.
The fact that our MLE contains a sum of logarithms of the observations should tell you that the given parametric distribution, if log-transformed, will yield a distribution for which the corresponding MLE will be a function of the (arithmetic) sample mean. In other words, let $Y = -\log X$ (we negate the log because when $0 < X < 1$, we want $Y > 0$). Then the density of $Y$ is $$f_Y(y) = f_X(e^{-y}) e^{-y} = (1+\theta)e^{-\theta y} e^{-y} = (1 + \theta) e^{-(1+\theta) y}, \quad y > 0.$$ This is an exponential distribution with rate $\lambda = 1+\theta$. It is easy to compute the MLE of $\lambda$ given a log-transformed sample $$y_i = -\log x_i, \quad i = 1, 2, \ldots, n:$$ the likelihood is $$\mathcal L(\lambda \mid y_1, \ldots, y_n) \propto \lambda^n e^{-\lambda \sum_{i=1}^n y_i} = \lambda^n e^{-n \lambda \bar y},$$ where $\bar y$ is the sample mean of log-transformed values. Then the log-likelihood is $$\ell(\lambda \mid y_1, \ldots, y_n) = n \log \lambda - n \bar y \lambda$$ which attains a maximum when $$\frac{d\ell}{d\lambda} \propto \frac{1}{\lambda} - \bar y = 0,$$ or $$\hat \lambda_{\text{MLE}} = \frac{1}{\bar y}.$$ Consequently, $$\hat \theta_{\text{MLE}} = -1 + \hat \lambda_{\text{MLE}} = -1 + \frac{1}{\bar y},$$ but we established that $\bar y$, being the sample mean of (negated) log-transformed values, is $$\bar y = -\frac{1}{n} \sum_{i=1}^n \log x_i.$$
To understand what is really going on, it helps to consider a concrete numeric example. Suppose we generate $n = 7$ observations from $X$ with $\theta = 3.14$. I have done this for you: $$(x_1, \ldots, x_7) = (0.86454, 0.839549, 0.890887, 0.886008, 0.867728, 0.958648, 0.8069).$$ Now, if we compute the MLE using the formula $$\hat \theta_{\text{MLE}} = -1 - \frac{n}{\sum_{i=1}^n \log x_i},$$ we first compute the negated log-transformed sample $$(y_1, \ldots, y_7) = (0.145558, 0.17489, 0.115538, 0.121029, 0.141877, 0.0422313, 0.214556).$$ Note that these observations are exponentially distributed with rate $\lambda = \theta + 1 = 4.14$. We take the sample mean: $$\bar y = 0.136526.$$ Then our MLE for $\theta$ is $$\hat \theta_{\text{MLE}} = -1 + \frac{1}{\bar y} = 6.32463.$$ Of course, the MLE is not very good because $n = 7$ is a small sample. If we do this for a large $n$, say $n = 10^5$, the performance is much improved: in Mathematica, we can perform such a simulation with the input
and the output I got was $3.14906$.