Particle model: first order reaction - Multiphase Reactor Modeling using PyMRM

Numerically solve the stationary spherically symmetric problem

\frac{1}{r^{2}}\frac{\partial}{\partial r}\left( r^{2}D\frac{\partial c}{\partial r} \right) - kc = \mathrm{div}\left( D\ \mathrm{grad}(c) \right) - kc = 0

with the boundary condition $c(R) = 1$ .

Questions:

Perform the spatial discretization (start with a uniform grid) and implement it in Python. Note that the spherical geometry is accounted for by the proper definition of the divergence operator.
Construct, implement, and solve the matrix-vector equation.
Compute the apparent reaction rate from the concentration gradient at the surface of the particle.
Does the result correspond to the analytical solution?
Investigate the effectiveness as a function of the Thiele modulus.
Consider the high Thiele modulus case ( $\phi\gg 1$ ) and improve the solution by using a spatial discretization that is refined near the wall.