# Question about coupled pde

**URL:** <https://fenicsproject.discourse.group/t/question-about-coupled-pde/3359>\
**Category:** Uncategorized\
**Created:** [May 29, 2020, 3:41am UTC](https://fenicsproject.discourse.group/t/question-about-coupled-pde/3359 "2020-05-29T03:41:42Z")\
**Posts on this page:** 1\
**Page:** 1

<div class="post-metadata">

**Author:** ![bbaek](https://avatars.discourse-cdn.com/v4/letter/b/e47c2d/32.png) [@bbaek](https://fenicsproject.discourse.group/u/bbaek)\
**Post date:** [May 29, 2020, 3:41am UTC](https://fenicsproject.discourse.group/t/question-about-coupled-pde/3359/1 "2020-05-29T03:41:42Z")

</div>

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.

![1](https://global.discourse-cdn.com/free1/uploads/fenicsproject1/original/2X/5/52c271be3075fd69810cc72df42b64fab5db4d99.png)

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.0_u1_exp(18\*(1-u2\*\*-1))  
#end def

u.interpolate(Constant((0.1, 1.0)))  
w = as\_vector([1.0])  
F1 =1/Pe_dot(grad(u\_1), grad(v\_1))dx - u\_1v\_1_ds(1) + v\_1_ds(1)   
+ dot(w, grad(u\_1))v\_1dx   
+ R(u\_1, u\_2)v\_1dx  
F2 =1/Bo_dot(grad(u\_2), grad(v\_2))_dx - u\_2_v\_2_ds(1) + v\_2_ds(1)  
+ dot(w, grad(u\_2))_v\_2_dx   
+ Hw\*R(u\_1, u\_2)_v\_2_dx  
F = F1 + F2  
solve(F == 0, u)