Problem while solving a minimal 1D Poisson-Nernst-Planck equation

Question

I am trying to solve the Poisson-Nernst-Planck equation in Python with the library FiPy. It is basically a set of equations that describes - in my example - the separation of two ion concentrations in a solution with a potential gradient. Like mixing salt in water and then applying a voltage-difference between both ends of the solution. The equations are:

And the boundary conditions:

Problem

I manage to get a converging solution, but it is not the solution I'd like. Specifically, I get that C_p and C_n always converge to be spatially constant, and is always linear. It seems like is indifferent to initial conditions. Even setting internal fixed point (based on this answer) turns to be linear with a break point, which doesn't help.

The solution I look for should be that is more like a sigmoid, and C_p and C_n get high values around the edges and low values in the center.

I really would appreciate help!

from fipy import *


# Constants
Lx = 10
nx = 1000
dx = Lx / nx # mm 

D = 1 
epsilon = 1

# FiPy
mesh = Grid1D(dx=dx, Lx=Lx)
x = mesh.cellCenters[0]

Cp = CellVariable(name="$C_p$", mesh=mesh, hasOld=True)
Cn = CellVariable(name="$C_n$", mesh=mesh, hasOld=True)
phi = CellVariable(name="$\phi$", mesh=mesh, hasOld=True)

viewer = Viewer((Cp, Cn, phi), limits={"ymax":5, "ymin":-5})

# Initial
Cn.setValue(0.5)
Cp.setValue(0.5)
# phi.setValue(6/(1+numerix.exp(-(x-4)))-3) # Sigmoid

# Boundry
Cp.faceGrad.constrain(0, mesh.facesRight)
Cp.faceGrad.constrain(0, mesh.facesLeft)
Cn.faceGrad.constrain(0, mesh.facesRight)
Cn.faceGrad.constrain(0, mesh.facesLeft)
phi.constrain(3, mesh.facesRight)
phi.constrain(-3, mesh.facesLeft)


# Equations
Cp_diff_eq = TransientTerm(coeff=1, var=Cp) == DiffusionTerm(coeff=D, var=Cp) + DiffusionTerm(coeff=D*Cp, var=phi)
Cn_diff_eq = TransientTerm(coeff=1, var=Cn) == DiffusionTerm(coeff=D, var=Cn) - DiffusionTerm(coeff=D*Cn, var=phi)
poission_eq = DiffusionTerm(coeff=epsilon, var=phi) == (ImplicitSourceTerm(var=Cn) - ImplicitSourceTerm(var=Cp))

equations = poission_eq & Cp_diff_eq & Cn_diff_eq

# Simulation
timestep = 0.01
time_final = 20
desired_residual = 1e-2

time = 0
i = 0
while time < time_final:
    phi.updateOld()
    Cp.updateOld()
    Cn.updateOld()

    residual = 1e10
    j=0    
    while residual > desired_residual:
        print(f"{i}-{j}")
        residual = equations.sweep(dt=timestep)
        j+=1
    
    if i % 10 == 0:
        viewer.plot()
    
    time_inc = timestep
    time += time_inc
    i += 1