I am looking at a few reaction-diffusion equations of the form
$\frac{dP}{dt} = D\left(\frac{d^2P}{dr^2} + \frac{2}{r}\frac{dP}{dr}\right) - a(P)$
I know the initial conditions and the boundary value at one end. I also know the function steadily falls and eventually hits zero at the second unknown boundary. To model this behaviour and indirectly find this unknown boundary, I coded a solver using an explicit finite difference scheme, with an added condition in the loop that changed forced any negative values to zero. This gave me the correct result, but due to the CFL condition I've had to use a tiny time step ($\Delta t = 0.001$) when I'm more interested in the system's behaviour over several hours or even days.
I was looking into using implicit finite difference methods or CN methods so I may increase the time step, but my (limited) understanding of this implies I'd have to solve a system of equations which would include the unknown boundary, which I will not know exactly. Is it possible to work around this, or will implicit methods fail? If it is possible to work around it, could anyone suggest a good method and how I would implement this? Thanks in advance, I'm quite new to numerical methods and would appreciate any suggestions.