Solving nonlinear equation with an embedded matrix equation

I want to find the solution to a nonlinear equation of the form u_{xx} + f(u,\lambda) + p= 0 where the unknowns are u(x) and the constant p. There is an addition volume constraint \int u(x)dx=1. For the leading order problem, f(u) is an explicit function of u and some parameters \lambda. After discretization using finite difference, I manage to solve it using packages like NLsolve.jl or even the full family of solutions (with different \lambda) using continuation in BifurcationKit.jl.

Suppose I want to tackle the higher order problem, now f_i=f(u_i) (the parameter \lambda is omitted for a sec just for simplicity) is no longer an explicit function of u (f_i, u_i are the discretized f(u) and u). To obtain f_i, I need to solve a matrix equation Af=b where the entries of the (sparse) matrix A and vector b depend on u. So the whole function I’m trying to solve looks something like this

function f(x, par)

u = @view u[1:end-1]
p = @view u[end]

uxx = finiteDiff(u, 2, par.order, par.dx)  # my custom function for F.D.

A, b = assemble_matrices(u, par)   # construct A and b using u(x) and given parameters
f = A \ b

residual = similar(x)
residual[1:end-1] = uxx .+ f .+ p
residual[end] = sum(u)/length(u) - 1

return residual

I’m aware of things like using non-allocating arrays. The above is just a high-level overview of the function f. Now the problem is that when I use the packages I mentioned above, they usually rely on automatic differentiation, and then very often I get StackOverFlow error inside ForwardDiff when it gets to the f=A\b line. Using finite difference for the Jacobian seems to work, but due to the matrix equation solve it is much slower my leading order problem. I just wonder if there are some general techniques for solving nonlinear equations like this? Do I just need to try my best to simplify the construction of matrix A, use as many non-allocating arrays as possible, or I just have to reformulate the problem to get rid of f=A\b?

Do you have an MWE that you can share? Nothing about A\b is bad for automatic differentiation, so what you have described is somewhat of a red herring. Whatever is giving the stackoverflow is something not in the description here. Some questions that would be helpful:

  • What are the types of A and b?
  • Can you isolate this A\b stackoverflow to be outside the nonlinear solver?