Code for https://stackoverflow.com/q/69267851/2019542

matlabish.py

from fipy import Grid1D, CellVariable,TransientTerm, DiffusionTerm, Viewer
from fipy.tools import numerix

nx = 100
Lx = 1.

mesh = Grid1D(nx=nx, Lx=Lx)
u1 = CellVariable(name = "solution variable u1", hasOld=True, mesh = mesh, value = 0.)
u2 = CellVariable(name = "solution variable u2", hasOld=True, mesh = mesh, value = 0.)

#initial conditions

u1.setValue(1)
#u2.setValue(0) #unnecessary

#boundary conditions

# u1.faceGrad.constrain([0.], where=mesh.facesLeft)  # no-flux is default BC
u2.constrain(0.0, where=mesh.facesLeft)  #Direction?
# u2.faceGrad.constrain([0.], where=mesh.facesLeft)  # no-flux is default BC
u1.constrain(1.0, where=mesh.facesRight) #Direction?


eqn0 = TransientTerm(var=u1) == DiffusionTerm(0.024, var=u1) - (numerix.exp(5.73*(u1-u2))-numerix.exp(-11.46*(u1-u2)))
eqn1 = TransientTerm(var=u2) == DiffusionTerm(0.170, var=u2) + (numerix.exp(5.73*(u1-u2))-numerix.exp(-11.46*(u1-u2)))
    
eqn = eqn0 & eqn1
#compute equations

viewer = Viewer(vars=(u1, u2), datamin=0., datamax=1.)
viewer.plot()

from steppyngstounes import PIDStepper

u1.updateOld()
u2.updateOld()

for step in PIDStepper(start=0., stop=2., size=0.0002):
    res = 1e100
    for sweep in range(1):
        res = eqn.sweep(dt=step.size)
    if step.succeeded(error=res / 0.01):
        u1.updateOld()
        u2.updateOld()
        viewer.plot()
        print(step.end, ",", res)
    else:
        u1.value = u1.old
        u2.value = u2.old