Numerical solution to parametrized second order ODE with nonuniform coefficients

I am trying to solve numerically the following second order linear ODE:

$$a \frac{\partial^2 u}{\partial x^2} + \frac{\partial u}{\partial x} \frac{\partial a}{\partial x} + b u =0$$,

on the domain $$[-L,L]$$ with boundary conditions $$u(-L)=u(L)=0$$

where $$a=\lambda c(x) + d(x)$$, and $$b = \lambda e(x) + f(x)$$.

$$c(x)$$, $$d(x)$$, $$e(x)$$, $$f(x)$$ are known functions of $$x$$ and $$\lambda$$ is a real parameter. It is a bit similar to Sturm-Liouville problems, but the parameter $$\lambda$$ is found in the coefficients of the derivatives. I don't know if any of the algorithms for that kind of problems apply, or if there is something else.

I am also interested if there is any method to determine $$\lambda$$ for which we have a nontrivial solution. Thank you!

You should be able to set this up as some sort of generalized eigenvalue problem like $$Au=\lambda Bu,$$ where $$Au = (du')'+fu$$ and $$Bu = -(cu')'-eu$$. I simplified the derivatives a bit since this formulation is more numerically stable. You should now be able to represent both as tridiagonal matrices using standard finite differences and any generalized eigenvalue problem solver will give you what you want.