Hi,
I am very new to Fenics and currently learning to use Fenics in my work. As a starting point, I am trying to replicate some literature. The problem is coupled PDE associated with simple chemical reactions.
I noticed that there are differences in solutions from the literature. It seems like the Robin boundary condition at x = 0 causing differences. I checked with other software too and was able to replicate the literature data. Can anybody know what seems to be a problem? Thanks in advance.

from fenics import *
import matplotlib.pyplot as plt

mesh = UnitIntervalMesh(10)
P1 = FiniteElement(‘P’, mesh.ufl_cell(), 2)
element = MixedElement([P1, P1])
V = FunctionSpace(mesh, element)
v_1, v_2 = TestFunctions(V)
u = Function(V)
u_1, u_2 = split(u)
color = MeshFunction(“size_t”, mesh, 0)
color.set_all(0)

bc1 = ‘x[0] == 0’
CompiledSubDomain(bc1, tol=1e-15).mark(color, 1)
ds = Measure(‘ds’, domain=mesh, subdomain_data=color)

Pe, Bo = Constant(10), Constant(10)
Hw = Constant(-0.05)

def R(u1, u2):
return 4.0u1exp(18*(1-u2**-1))
#end def

u.interpolate(Constant((0.1, 1.0)))
w = as_vector([1.0])