# Implementation of Neumann boundary condition with method of lines - 1D diffusion/reaction equation

I am solving the monodimensional diffusion/reaction equation by discretization using the well-known method of lines

${\partial c\over\partial t}=D{\partial^2 c\over\partial x^2}+ r\text{ for t>0 and -L<x<L} \quad[1]$
$c(x,0)=1 \quad[2], \quad \text{ initial condition}$
${dc\over dx} \mid_{(0,t)} =0\ \text{ for t>0,} \:[3], \text{ symmetry at the center}$
$D{dc\over dx} \mid_{(L,t)} =k_L(C_L-C\mid_{(L,t)}) \quad[4]$

This last equation represents that the diffusion at the boundary equals the external mass transfer but I am having trouble to transform it to the method of lines. The discretization of the system becomes:
${dc_i\over dt} =D{2 c_{i+1}-c_i\over \Delta x^2}+r \text{ for i=1} \quad[5]$
${dc_i\over dt} =D{ c_{i+1}-2c_i+c_{i-1}\over \Delta x^2} +r\text{ for i = 2...N-1} \quad[6]$

And then, for the last condition I am not sure how to implement it. Discretizing equation [4], I obtain $C_i ={k_L\:c_L + {D\over \Delta x} c_{i-1} \over k_L + {D\over \Delta x} } \text{ for i = N} \quad[7]$

I don't think that this is correct given that $C_N$ does not depend on time, is calculated indirectly at each time step from $C_{N-1}$. And that I cannot obtain the derivative for ${dc_{N-1}\over dt}$. But I cannot see how that condition should be discretized. Thank you for your help.

• Equation [7] looks fine to me. – David Ketcheson Dec 29 '15 at 22:20
• Then you would first solve eq. [7] and then eq. [6] to obtain the derivative and N-1? – Toulousain Dec 29 '15 at 23:21
• [6] is continuous in time, so you still have some discretizing to do. You can view the full discretization as a coupled system or use [7] to eliminate an unknown. – David Ketcheson Dec 30 '15 at 2:37
• Thank you very much, Peter, for your answer. Your answer helped me in brushing up my knowledge. I would like to mention one point. When I applied the 2nd Boundary Condition, I acheived this: dc(N)/dt = D[-c(N)+c(N-1)]/dx^2 + kL/dx*[CL-c(N)] + r – Iftekhar Aug 3 '17 at 8:24

Although in principle your method can work, you might need a quite fine mesh to obtain a good approximation, because it is like you solve it by MOL on a domain $(h,L-h)$ and with only a first order accurate scheme.
Let the index $i$ denote in which space coordinate $x_i$ you approximate your equation. I prefer $x_i = i h$ for $h=L/N$ and $i=0,1,\ldots,N$. Your aim should be to construct a system of ODEs for unknown functions $c_i(t) \approx c(x_i,t)$ for $i=0,1,\ldots,N$, it means including the concentrations at boundary nodes.
Using central difference to approximate $dc/dx$ in Neumann boundary conditions and using it to eliminate auxiliary functions $c_{-1}(t)$ and $c_{N+1}(t)$ you should obtain $$\frac{dc_0}{dt}=D \frac{-2c_0+2c_1}{dx^2}+r \quad [5]$$ and $$\frac{dc_N}{dt}=D \frac{-2c_N+2c_{N-1}}{dx^2}+\frac{2k_L}{dx}(C_L-c_N)+r \quad [7]$$