Solver for Elliptic PDE with Dirichlet and Neumann BCs

Hello,

I am interested in developing a finite difference solver for a 2D elliptic PDE (~ Darcy flow) from (https://arxiv.org/pdf/2205.09322.pdf). Do you have recommended packages to solve this PDE?

The PDE is given by:

-\text{div}(e^{u(x)} \nabla u(x)) = f(x), \text{for } x=(x_1, x_2) \in [0,1] \times [0,1]
with boundary conditions:
v(x_1, 0) = 100, \frac{\partial v}{\partial x_1}(1, x_2) = 0, -e^{u(x)} \frac{\partial v}{\partial x_1}(0, x_2) = 500, \frac{\partial v}{\partial x_2}(x_1, 1) = 0,
and forcing term
f(x) = f(x_1, x_2) = \begin{cases} 0 & 0 \leq x_2 \leq \frac{4}{6},\\ 137 & \frac{4}{6} < x_2 \leq \frac{5}{6},\\ 274 & \frac{5}{6} < x_2 \leq 1 \end{cases}

We assume that the log diffusion term is given by
u(x_1, x_2) = \sum_{i=0}^{19} \sum_{j=0}^{19} u_{i,j} \varphi_{i,j}(x_1, x_2),

with \varphi_{i, j}(x_1, x_2) = \cos(i \pi x_1) \cos(j \pi x_2).

A coarse mesh ~ 20 \times 20 grid points is sufficient =)

I would appreciate your help.