Hi, I need to model a reaction of rate constant k and find the concentration of both reactant Cr and product Cp that follow standard diffusion-reaction equations i.e.:

d(Cr)/dt = d^2(Cr)/dx^2 - k(Cr)

d(Cp)/dt = d^2(Cp)/dx^2 + k(Cr)

for just the first equation, I have successfully implemented the Crank-Nicolson method (euler method is not accurate enough) but I run into difficulties when I try to do the second one - I don't know how to deal with having both Cr and Cp... is this even possible?

Any help would be greatly appreciated.