Hi, I’m currently solving the diffusion problem of a plate, but when I apply a boundary on both sides, I always get a non-convergent result，I hope you can give me some suggestions, here is an example of my code：

```
from fenics import *
import numpy as np
import matplotlib.pyplot as plt
# 定义已知表达式
kB = Constant(8.6e-5)
r = Constant(0.3)
Sext = Constant(0.0)
fai = Constant(1e24)
f = 0.002
u_0 = Expression('0', degree=1)
t_total = 10
num_steps = 500
dt = t_total / num_steps
def D(T):
return 4.17e-7*exp(-0.39/(kB*T))
def KrW(T):
return 3.2e-15*exp(-1.16/(kB*T))
def KrCuCrZr(T):
return 2.9e-14*exp(-1.92/(kB*T))
def S(T):
return 9e-3*exp(-1.04/(kB*T))
mesh = Mesh('10_10_100s_slab/mesh.xml.gz')
W = FunctionSpace(mesh, 'P', 1)
V = FunctionSpace(mesh, 'P', 1)
class BoundaryX0(SubDomain):
def inside(self, x, on_boundary):
return on_boundary and abs(x[0] - 0.00) < DOLFIN_EPS
class BoundaryX1(SubDomain):
def inside(self, x, on_boundary):
return on_boundary and abs(x[0] - 0.01) < DOLFIN_EPS
class BoundaryY0(SubDomain):
def inside(self, x, on_boundary):
return on_boundary and abs(x[1] - 0.00) < DOLFIN_EPS
class BoundaryY1(SubDomain):
def inside(self, x, on_boundary):
return on_boundary and abs(x[1] - 0.01) < DOLFIN_EPS
bx0 = BoundaryX0()
bx1 = BoundaryX1()
by0 = BoundaryY0()
by1 = BoundaryY1()
boundary_markers = MeshFunction('size_t', mesh, mesh.topology().dim()-1) # 创建网格函数
#boundary_markers.set_all(9999)
bx0.mark(boundary_markers, 0)
bx1.mark(boundary_markers, 1)
by0.mark(boundary_markers, 2)
by1.mark(boundary_markers, 3)
ds = Measure('ds', domain=mesh, subdomain_data=boundary_markers)
T = Function(W)
c1 = Function(V)
c_n = project(u_0, V)
v = TestFunction(V)
K1 = KrW(T)*S(T)**2
F = c1*v*dx + dt*D(T)*dot(grad(c1),grad(v))*dx - c_n*v*dx - dt*Sext*v*dx - dt*((1-r)*fai-KrW(T)*c1**2)*v*ds(3) + dt*(KrCuCrZr(T)*c1**2)*v*ds(2) - dt*(KrW(T)*c1**2-K1*f)*v*ds(0)
timeseries_T = TimeSeries('10_10_100s_slab/temperature_series')
progress = Progress('Time-stepping', num_steps)
t = 0
bcs=[]
J = derivative(F, c1)
problem = NonlinearVariationalProblem(F, c1, bcs, J)
solver = NonlinearVariationalSolver(problem)
prm = solver.parameters
prm['nonlinear_solver'] = 'newton'
prm['newton_solver']['absolute_tolerance'] = 1e-8
prm['newton_solver']['relative_tolerance'] = 1e-5
prm['newton_solver']['maximum_iterations'] = 100
prm['newton_solver']['relaxation_parameter'] = 1.0
plt.figure(figsize=(12,8))
for n in range(num_steps):
t += dt
timeseries_T.retrieve(T.vector(), t)
solver.solve()
plot(c1)
c_n.assign(c1)
set_log_level(LogLevel.PROGRESS)
progress += 1
#set_log_level(LogLevel.ERROR)
```

Be grateful！