I'm trying to use a finite difference method to solve the steady-state neutron diffusion equation in a nuclear reactor:
$$ D(x) \nabla^2 \phi(x) + \left( \frac{\nu(x)}{k} \Sigma_f(x) - \Sigma_a(x) \right) \phi(x) = 0 $$
When the material properies $D$, $\nu$, $\Sigma_f$ and $\Sigma_a$ are uniform through the reactor, this reduces to
$$ \nabla^2 \phi(x) + B^2 \phi(x) = 0 $$
where $B^2 = \frac{\frac{\nu}{k} \Sigma_f - \Sigma_a}{D}$. I can create the matrix equation
$$ A \mathbf{u} = \lambda \mathbf{u} $$
and the solution I'm looking for is the eigenvector corresponding to the smallest magnitude eigenvalue $\lambda$. That eigenvalue is the $-B^2$ required for the solution to be valid, and using $k=1$ (as is true for steady-state conditions), the required material properties of the reactor can be determined.
However, I don't know what matrix equation to construct when the material properties vary. In practice I want to have a number of regions each with different properties. When solving analytically for a one-dimensional reactor with two regions, the conditions at the interface are $\phi_1(x_{int}) = \phi_2(x_{int})$ and $D_1 \frac{d\phi_1}{dx}(x_{int}) = D_2 \frac{d\phi_2}{dx}(x_{int})$. Ideally, the properties of all but one of the regions would be defined, with the properties of the remaining region being determined by the conditions required to make $k$ equal to 1, as in the uniform properties case.
Please can someone point me in the right direction? Thanks very much for any help.