Let $X_1,\ldots,X_n$ be a sample of size $n$, all with pdf $f(x)=\frac{x+1}{\theta(\theta+1)}e^{\frac{-x}{\theta}}\Bbb I_{(0, \infty)}(x)$ where $\theta > 0$.
We will drop the indicator function to make the work a little shorter, although we acknowledge that $x>0$ for future calculations.
I am asked to find the MLE of $\theta$.
So $$\mathscr L(\theta, x_1,\ldots,x_n)=\prod_{i=1}^nf(x_i)=\prod_{i=1}^n\frac{x_i+1}{\theta(\theta+1)}e^{\frac{-x_i}{\theta}}=\frac{\prod_{i=1}^n(x_i+1)}{\theta^n(\theta+1)^n}e^{\frac{-\sum_{i=1}^n x_i}{\theta}}$$
We now take the log-likelihood, and we get
$$\ln(\mathscr L)=\ln \left(\frac{\prod_{i=1}^n(x_i+1)}{\theta^n(\theta+1)^n}e^{\frac{-\sum_{i=1}^n x_i}{\theta}}\right) = \ln\left(\prod_{i=1}^n (x_i+1)\right) - n\ln(\theta)-n\ln(\theta+1) - \frac{\sum_{i=1}^n x_i}{\theta}=\ln \left( \prod_{i=1}^n(x_i+1)\right)-n\ln(\theta) - n\ln(\theta+1)-\frac{n\overline X}{\theta}$$
Noting that $\sum_{i=1}^n x_i=n\overline X$, where $\overline X$ is the sample mean.
Now we take the derivative of $\ln(\mathscr L)$ with respect to $\theta$ and let it equal $0$.
$$\frac{\partial}{\partial \theta}\ln(\mathscr L)=-\frac{n}{\theta}-\frac{n}{\theta+1}+\frac{n\overline X}{\theta^2}=0$$
$$\frac{n\overline X}{\theta^2}=\frac{n}{\theta}+\frac{n}{\theta+1}$$
$$\frac{n\overline X}{\theta^2}=\frac{2n\theta+n}{\theta(\theta+1)}$$
$$\frac{\overline X}{\theta}=\frac{2\theta+1}{\theta+1}$$
$$\frac{\overline X}{\theta}=\frac{2\theta+2}{\theta+1}-\frac{1}{\theta+1}$$
$$\frac{\overline X}{\theta}=2-\frac{1}{\theta+1}$$
$${\overline X}=2\theta-\frac{\theta}{\theta+1}$$
$$(\theta+1){\overline X}=2\theta(\theta+1)-\theta$$
$$(\theta+1){\overline X}=2\theta^2+2\theta-\theta$$
$$(\theta+1){\overline X}=2\theta^2+\theta$$
$$-(\theta+1){\overline X}+2\theta^2+\theta=0$$
$$2\theta^2+\theta(1-\overline X)-\overline X=0$$
At this point I tried using the quadratic formula but it doesn't make much sense how it works out, then the MLE ends up being
$$\frac{(\overline X-1)\pm \sqrt{(1-\overline X)^2+8\overline X}}{4}$$
I have a feeling this doesn't make much sense. Where did I go wrong?
Also, I would then need to find the MLE of $\theta(2\theta+1)/(\theta+1)$, but I have a feeling it simply follows by plugging in the MLE found for $\theta$ and then simplifying.
Your last sentence follows from the "invariance property" of MLEs. (It probably ought to be called equivariance rather than invariance.)
This family of distributions is somewhat weird: It's a weighted average of two gamma distributions with different shape parameters and the same scale parameter, with the two weights depending on the scale parameter. Maybe that accounts for any seeming weirdness in the answer.
It should be remembered that an MLE is not defined simply as any solution of the likelihood equation. It's a global maximum point. The value of a function's derivative can be $0$ at points other than global maxima.
Going through your work somewhat quickly, I don't see anything wrong in it, except that you're a bit overly credulous if you think it shows that both solutions of the quadratic equation are MLEs. MLEs (unlike unbiased estimators and method-of-moments estimators) are never outside of the parameter space (or its closure, if you don't insist of the parameter space being a closed set (and if you don't, then you might say the MLE doesn't exist if it's on the boundary)).
I once attended a talk by Herbert Robbins in which he said that after much work he finally concluded that a number he sought in a problem of statistics satisfied a certain quadratic equation, and then it took him eight month of further research to decide which of the two solutions it was.