Does wing with bell lift distribution has negative lift at tips?
1. Summary
As your first plot shows, the optimal Bell-Shaped Lift Distribution (BSLD) is positive throughout its span for an untwisted planform. What's different is that the downwash near the wingtip is negative for the BSLD, thereby giving rise to an apparent thrust and decreasing the overall induced drag compared to the Elliptical-Shaped Lift Distribution (ESLD).
However, for a non-optimal BSLD, you may get negative lift near the wingtips. And when there is geometric/aerodynamic twist, then you will almost certainly get negative lift over some portion of the span at AOAs lower than your design point.
2. Reasoning
Inspired by an Aerocrafty blog (be cautious that the derivation in that article contains some minor errors), we can derive an approximate BSLD without resorting to complex analysis. From Lifting-Line Theory, section lift ($L^{'}$) and upwash ($w$) are given as follows:
$$L^{'}(y_0) = \rho_\infty V_\infty \Gamma(y_0)$$
$$w(y_0) = -\frac{1}{4\pi} \int_{-b/2}^{b/2}\frac{(d\Gamma/dy)}{y_0-y}dy$$
where $y_0$ is the lateral span coordinate, $b/2$ is half span, $\Gamma(y_0)$ is the circulation distribution of the bound vortex, $\rho_\infty$ is the free-stream density, and $V_\infty$ is the free-stream airspeed.
Through a transformation of coordinates (from Cartesian to polar):
$$y=-\frac{b}{2}\cos{\theta}$$
We can write any general circulation distribution as a Fourier sine series:
$$\Gamma(y) = \Gamma(\theta) = \sum_{n=1}^{\infty} {A_n\sin{n\theta}}$$

For an ESLD, the circulation distribution is only the first term, with $A_1=\Gamma_{e_0}$, where $\Gamma_{e_0}$ is the circulation at wing root:

$$\Gamma_e(\theta) = \Gamma_{e_0}\sin{\theta}$$

For a BSLD, we will only keep the first two terms and truncate the rest (noting that due to symmetry, only odd terms remain):

$$\Gamma_b(\theta) \approx A_1\sin{\theta} + A_3\sin{3\theta}$$
Now we can solve for the unknown coefficients $A_1$ and $A_3$ in the BSLD distribution by requiring that a comparative ESLD must have the same total lift ($L$) and root bending moment ($B$):
$$L = \int_{-b/2}^{b/2} {L^{'}dy} = L_e = L_b \tag{a}$$
$$B = \int_{-b/2}^{b/2} {yL^{'}(y)dy} = B_e = B_b \tag{b}$$
That's two equations and two unknowns. The solutions are:
$$A_1 = \frac{b_e}{b_b}\Gamma_{e_0}$$
$$A_3 = \frac{5\Gamma_{e_0}}{3} \frac{b_e}{b_b} \left( \frac{b_e}{b_b}-1 \right)$$
where $b_e$ is the wing span in the ESLD case satisfying (a) and (b), while $b_b$ is the wing span in the BSLD case. That is, we are expressing BSLD as a ratio of span ratio compared to an Elliptical case with the same total lift and root bending moment. Notice that when $b_e/b_b = 1$, we retrieve the ESLD exactly.
2.1 Validation
In the beginning, I mentioned that this derivation is an approximation. But it's not a bad one. Calculating the total induced drag with this formulation, we get:
$$D_i = \frac{2L^2}{\rho_\infty \pi b^2 V_\infty^2} \left[\frac{25}{3} \left( \frac{b_e}{b_b} \right)^4 - \frac{50}{3} \left( \frac{b_e}{b_b} \right)^3 + \frac{28}{3} \left( \frac{b_e}{b_b} \right)^2 \right]$$
Compare this with the exact result from NACA TN-2249:
$$D_i = \frac{2L^2}{\rho_\infty \pi b^2 V_\infty^2} \left[8\left( \frac{b_e}{b_b} \right)^4 - 16\left( \frac{b_e}{b_b} \right)^3 + 9\left( \frac{b_e}{b_b} \right)^2 \right]$$
2.2 Results
Time to get back to the question we began with, can the lift (or circulation) distribution be negative on any part of the span? The plot below shows the circulation distribution of our approximate solution for various ratios:

For our approximate formulation, the most optimal BSLD occurs at $b_e/b_b = 0.8$ ($b_b/b_e = 1.25$ in the plot). As you can see, it's all positive up to the wing tip.
However, as we keep increasing the span (while maintaining total lift and root bending moment), at some point the distribution goes negative near the wingtip. Eventually the bending moment at certain points along the span will exceed that at the wing root (Ref NACA TN-2249).