To complement Soner's answer, I want to clarify why we cannot naively simply plug in $\dot{\theta}=l/mr^2$ in the Lagrangian. This is discussed in a different context in exercise 1.5 of Henneaux, Marc, and Claudio Teitelboim. Quantization of Gauge Systems. Princeton University Press, 1992.
First of all, let us recall that if we want to obtain the equations of motion for an orbit starting at $(r_i,\theta_i)$ and finishing at $(r_f,\theta_f)$, we want to minimize the action
$$S(r,\theta)=\int_{t_i}^{t_f}\text{d}t\left(\frac{1}{2}m\dot{r}^2+\frac{1}{2}mr^2\dot{\theta}^2-U(r)\right)$$
on the space $\mathcal{E}$ of trajectories $r=r(t)$ and $\theta=\theta(t)$ satisfying $r(t_i)=r_i$, $\theta(t_i)=\theta_i$, $r(t_f)=r_f$, and $\theta(t_f)=\theta_f$. What one can show, say from Noether's theorem, is that any minima of this action satisfies that
$$l=mr^2\dot{\theta}$$
is a constant. This would imply that at the minima of the action, $\dot{\theta}=l/mr^2$ for some constant $l$. By integrating both sides, we see that this constant, by virtue of the boundary conditions of the space we are working with, is, in fact, a functional of the trajectory $r=r(t)$
$$l(r)=\frac{m\Delta\theta}{\int_{t_0}^{t_f}\text{d}t\,\frac{1}{r^2}},\quad\Delta\theta=\theta_f-\theta_i.$$
You see, the difficulty here is one of language. While $l$ is indeed a constant over time, it is not a constant on the space of trajectories $\mathcal{E}$.
We can now proceed keeping this in mind. On the minima, we know that
$$\int_{t_0}^{t_f}\text{d}t\,\frac{1}{2}mr^2\dot{\theta}^2=\int_{t_0}^{t_f}\text{d}t\,\frac{l(r)^2}{2mr^2}=\frac{l(r)^2}{2m}\int_{t_0}^{t_f}\text{d}t\,\frac{1}{r^2}=\frac{\Delta\theta}{2}l(r).$$
Therefore, the minima is also obtained by minimizing on $X$ the function
$$\tilde{S}(r,\theta)=\int_{t_i}^{t_f}\text{d}t\left(\frac{1}{2}m\dot{r}^2+\frac{l(r)^2}{2mr^2}-U(r)\right)=\int_{t_i}^{t_f}\text{d}t\left(\frac{1}{2}m\dot{r}^2-U(r)\right)+\frac{\Delta\theta}{2}l(r)$$
This is a well-posed problem, although it is not one that can be solved using Euler-Lagrange equations. The problem is that this action is not local, since the functional $l(r)$ is not the integral of some density. One can however proceed noting that
$$\delta l=\frac{2m\Delta\theta}{\left(\int_{t_i}^{t_f}\text{d}t\,\frac{1}{r^2}\right)^2}\int_{t_i}^{t_f}\text{d}t\,\frac{1}{r^3}\delta r=\frac{2l(r)^2}{m\Delta\theta}\int_{t_i}^{t_f}\text{d}t\,\frac{1}{r^3}\delta r.$$
We conclude
$$\delta \tilde{S}=\int_{t_i}^{t_f}\text{d}t\,\left(-m\ddot{r}-U'(r)+\frac{l(r)^2}{mr^3}\right)\delta r,$$
which gives the correct equations of motion
$$m\ddot{r}=-U'(r)+\frac{l^2}{mr^3}.$$
Now, while this clarifies the subtleties of replacing $\dot{\theta}=l/mr^2$ in the action, it doesn't bring any new light on how one can obtain from the action $\tilde{S}$ the equivalent action
$$S(r)=\int_{t_i}^{t_f}\text{d}t\,\left(\frac{1}{2}m\dot{r}^2-\frac{l^2}{2mr^2}-U(r)\right),$$
for a truly constant $l$, of course, other than noting that they both give the same minima. There is another procedure, shedding light into the role of boundary conditions, which does give the correct action.
Since adding boundary terms preserves the equations of motion, we can instead consider the action
$$\bar{S}(r,\theta)=\int_{t_i}^{t_f}\text{d}t\left(\frac{1}{2}m\dot{r}^2+\frac{1}{2}mr^2\dot{\theta}^2-U(r)\right)-[mr^2\theta\dot{\theta}]_{t_i}^{t_f}\,=\int_{t_i}^{t_f}\text{d}t\left(\frac{1}{2}m\dot{r}^2+\frac{1}{2}mr^2\dot{\theta}^2-U(r)-\frac{\text{d}}{\text{d}t}(mr^2\theta\dot{\theta})\right).$$
Indeed, the variation of such an action is of the form
$$\delta\bar{S}=\int_{t_0}^{t_f}\text{d}t\,\left((\dots)\delta r+(\dots)\delta\theta\right)+\left[(\dots)\delta r+(\dots)\delta{\dot{\theta}}\right]_{t_0}^{t_f}.$$
Thus, the variational problem for $\bar{S}$ is only well-posed if we consider spaces of trajectories where we fix boundary conditions for $\dot{\theta}$ instead of $\theta$. In such spaces, $l$ is truly a constant, e.g.
$$l=mr_i^2\dot{\theta}_i,$$
and we can consider the equivalent action in which $\dot{\theta}=l/mr^2$
$$\bar{S}(r,\theta)=\int_{t_i}^{t_f}\text{d}t\left(\frac{1}{2}m\dot{r}^2+\frac{l^2}{2mr^2}-U(r)\right)-l[\theta]_{t_i}^{t_f}\,=\int_{t_i}^{t_f}\text{d}t\left(\frac{1}{2}m\dot{r}^2+\frac{l^2}{2mr^2}-U(r)-l\dot{\theta}\right)=\int_{t_i}^{t_f}\text{d}t\left(\frac{1}{2}m\dot{r}^2-\frac{l^2}{2mr^2}-U(r)\right)=S(r).$$
In particular, we now obtain the correct effective potential term.
