Understanding the Born-Mayer binding potential for ionic crystals what was supposed to be a simple question, turned out to be a conundrum. I am asked to plot the Born-Mayer potential energy for a single pair of positive and negative ions. The potential energy is written thusly 
$$U(r)=\lambda e^{-r/\rho} - \frac{e_{0}^{2}}{r} $$
That $-1/r$ dependence is due to the colombic attraction between the positive and negative ion, and that exp(-r) repulsion comes from the Pauli Exclusion Principle. If I take the negative gradient of each, it checks out, since one results in a vector pointing towards the origin and the other away from it.
When I plot $U(r)$ however, the resulting graph looks nothing like a binding potential. That $-1/r$ dominates the expression for small enough $r$. The graph presents a local maximum instead of a local minimum.
If I flip the signs, everything seems to work out. Except that now I have changed the roles of each term, colombic attraction is now repulsion.
This is the graph I made: https://www.desmos.com/calculator/pmjjkyvept
I must be misunderstanding something.
 A: Your original starting point was correct. The Coulomb potential is attractive, while Born-Mayer potential is repulsive (positive $\lambda$). It is also correct that the sum of the two terms results in a short range diverging attraction at $r=0$, even if for some values may display a minimum (see for instance the values $\lambda=9$, $\rho=0.5$, with $ e_0=1$).
The key point (not always stressed enough when introducing the Born-Myer term) is that the potential is intended to be used *only in the region of the minimum and at large distencs. Absolutely not close to the maximum or on the left of it. Which is not a problem for minimization with respect to the nearest neighbor distances or for evaluating the dynamical matrix in the case of vibrational spectra. But, for example, when one has to simulate liquid alkali halides, some additional repulsion at short distances is introduced to avoid the possibility of a "collapse" of the system.
A: The constants are out of scale. The term $e_o$ for example has to be framed to consider the charge on an electron times the oxidation states of the two ions.
To compensate the mistake, use this equation in its more general form.
$$ U(r) = A \exp\left(-r/r_o\right) - B \frac{e_o^2}{r} $$
This will allow you to turn off the repulsive or attractive potentials by setting $A = 0$ or $B = 0$ respectively.
An potential curve illustrative of the shape is obtained with $A = 80$, $B = 10$, $r_o = 0.5$, and $e_o = 1$. This is not valid by physical scale to any real system. It does however show the minimum in the potential well, the hard-wall repulsion as $r$ is less than this value, and the extension to $U \rightarrow 0$ as $r \rightarrow \infty$ going from $U < 0$.
A better approach for any potential equation is to substitute the physical parameters as they should be. Alternatively, when that is not possible, find the potential $U_o$ at $r = r_o$ for the given equation and normalize to plot $U/U_o$. By example then, try this equation with $\Phi = U/U_o$ and $\rho = r/r_o$.
$$ \Phi = \left(\frac{1}{(1/b) - 1}\right)\left(\frac{1}{\rho} - \frac{1}{b} \exp[b (1 - \rho)]\right) $$
With $r_o = 1$ and $b = 10$, the plot of $\Phi$ (blue) and the plot of $F/F_{max}$ (green) are shown below.

A: I've read Jeffrey's and Giorgio's answers, thank you for your help. I would like to give something in return. Specifically, this how I determined an expression for $U/U_0$ for this particular potential. I will use the potential Jeffery stated, which is modified from the question to account for the oxidation of the ions.
$$U(r) = A\exp(-r/l) - B\frac{ e_0^{2}}{r}\tag1$$
We start by imposing that $U'(r_0)=0$. By doing so, we obtain a transcendental equation. Since we cannot obtain a closed-form for $r_0$, we'll need to get creative. We can manipulate the equation to obtain either one of these two results:
$$A\exp(-r_0/l) = lB(\frac{e_0}{r_0})^2 \tag2$$
$$B\frac{e_0^2}{r_0} = \frac{r_0}{l}A\exp(-r_0/l)\tag3$$
Note that they closely resemble the terms being summed in $U(r)$. We define $U_0 = U(r_0)$ and by substitution, we can obtain two expressions for $U_0$.
$$U_0 = B\frac{e_0^2}{r_0}(\frac{1}{\rho_0}-1) \tag4$$
$$U_0 = A\exp(-\rho_0)\left( 1-\rho_0 \right) \tag5$$
Here, I have defined $\rho_0 = r_0/l$. By hindsight, we know that there are two local extrema, one being a local maximum and the other a local minimum. As stated by GiorgioP, the model is relevant from the neighboorhood of the local minimum to large distances. As such, we now impose the condition that $U''(r_0)>0$, more explicitly
$$\frac{A}{l^2}\exp(-r_0/l) - 2B\frac{e_0^2}{r_0^3}>0 \tag6$$
By substituting equation (3) into (6), it becomes a simple matter of simplifying the inequation to obtain $\rho_0>2$. This implies that $U_0$ is negative. Finally, we divide the potential energy by $U_0$. The first term in $U(r)$ is divided by (5) and the second by (4). Define $\rho = r/l$, make all the relevant substitutions and the result should resemble the following equation
$$U(r)/U_0 = -\frac{1}{1-\rho_0}\left( e^{-(\rho-\rho_0)}-\frac{\rho_0^2}{\rho} \right) $$
The resulting graph can be seen here: https://www.desmos.com/calculator/grzz7bd9qs
