I am trying to solve a set of coupled, nonlinear ODEs. The only dependent variable is a 1-dimensional spatial coordinate, let's call it $x$. For now, I've managed to approximate away some of the coupling between the equations, such that the first one doesn't depend on the second one at all. I solved it with FEM without major difficulty: boundary conditions met, residuals look good, physically meaningful result.
So now I'm left with one nonlinear, second-order ODE where the nonlinear coefficient functions are completely determined by the FEM solution, something like
$$ \begin{cases} A(x) + B(x) f(x) + C(x) f'(x) + D(x) f''(x) = 0 \\ f(0) = f(1) = 0 \end{cases} $$
where $A$, $B$, $C$, and $D$ are known, at least at the resolution of the FEM mesh.
The problem is that, using the same basic code as before, the solution appears to violate the boundary conditions entirely, blowing up as $x \rightarrow 0$ instead of approaching zero (or even any finite value, it doesn't seem to matter what I choose). Thinking this was a problem with my code, I experimented with a finer mesh and with finite differences instead of finite elements, with no success. It finds the same nonsensical result every time.
Finally, I put my numerical data for the coefficients into Mathematica and attempted to solve the above BVP. I let it run for quite a while (maybe 15-30 minutes) before giving up. It never found a solution.
So finally, my question: Is it possible that this equation simply can't be solved with these boundary conditions? For example, I tried changing $f(1) = 0 \rightarrow f'(0) = 0$ and the solution was found almost instantaneously by Mathematica, but it still violated the new boundary condition; somehow, that didn't trigger an error, but it made me suspicious that something "bigger" might be going on here.
I'm looking for some mathematical intuition behind why the numerical methods might fail in these cases and how to mitigate these problems.
Update
This equation is linear (correcting here so that comments still make sense). I think that's not the root of my problem, though; I tried using a linear solver with the proper stiffness matrix and load vector to solve $-f''(x) = g(x)$, which should be EASY; here, $g(x)$ is the FEM solution to the first equation. The solution meets boundary conditions but is "jagged" (think sawtooth patterns overlaid on a smooth function). Now I think perhaps there's something wrong with substituting the FEM solution for the first equation into the second equation, but I'm not sure why.
Not giving up yet.
Update 2
I think I can safely rule out a problem with this code; I am generating the source files with SymPy, so there shouldn't be any human errors associated with changing functional forms, etc. of the equations I want to solve. As an example of what I'm seeing now, see figure below. The only intuition I have for this is that the code "tries" to meet boundary conditions (which are homogeneous Dirichlet), but the solution, for whatever reason, fundamentally wants to blow up to $-\infty$. So we end up with this "compromise" solution that straddles the two possibilities.
Is this even remotely correct? It doesn't seem to matter what I change at this point, and refining the mesh just increases the frequency of the sawtooth form.