I'm trying to solve the Poisson-Boltzmann equation with the Finite element approach:
$\begin{matrix} \frac{\partial^2 y}{\partial x^2}=\sinh{y}, & \\ \frac{\partial y}{\partial x} = a,b & at ~ x=0,h \end{matrix}.$
So I write:
$\int_{x_{i}}^{x_{i+1}} \left(\frac{\partial^2 y}{\partial x^2}-\sinh{y}\right) \omega (x) dx$
with the weight function as: $ \tilde{y}=\omega_{1} y_1 +\omega_{2} y_2$ and $\omega_1=\frac{x_{i+1}-x}{x_{i+1}-x_{i}}$, $\omega_2=\frac{x-x_{i}}{x_{i+1}-x_{i}}$.
and I get a system of equation :
$K_{i,j}y_{j}=F_{i}$
where:
$K_{i,j}=\int_{x_{i}}^{x_{i+1}} \frac{\partial y}{\partial x} \frac{\partial \omega}{\partial x} dx$
and ,
$F=\frac{\partial y}{\partial x} \omega \mid_{x_{i}}^{x_{i+1}}+\int_{x_{i}}^{x_{i+1}}\sinh{y}\omega (x) dx$
Finally using the Newton-Raphson method to deal with the non-linearity : $y^{it+1}=y^{it} -[K_{T}]^{-1} (K_{i,j} y^{it}-F(y^{it})$
where $K_T=K_{i,j}-\frac{\partial F_{i}}{\partial Y_{j}}$ and '$it$' is for the iteration.
However the results that I obtain do not match the analytical solution and I'm not sure what is the problem, is it because of the linear element and should I switch to cubic Hermite trial function or I'm missing something else?
comparison of results and analytical solution
Edit: I had modified the constant matrix according to this: link by writing $F=F-mean(F)$, as I thought imposing a pure Neumann boundary condition requires this modification. But I think this is not the case for Poisson-Boltzmann equation, since after removing it, I get acceptable accordance with the analytical solution. new results.
Now a new problem emerges, which is by lowering the domain size ($h$) below a value (say, 0.001), the $K_{T}$ matrix condition number gets high (~ $10^{16}$) and the numerical system breaks down. Is there a modification I can implement to circumvent this issue?