I have a reaction-diffusion equation in 2 dimensions of the typical form:

$\frac{\partial u}{\partial t} = D\nabla^2u - \Phi(u(x))$

I want to stress that $\Phi(u(x))$, is *not* a constant, but depends on the location, where it can either be a source or a sink term for the field $u$. The rate of consumption for the sink term here depends on the level of the field $u(x)$. The source of the field is constant.

At the moment, I have been assuming that the field is in steady-state at each time step, so have been solving the nonhomogeneous Poisson equation of the form:

$D\nabla^2u = \Phi(u(x))$

I have been using a finite difference approximation (where alterations are made to the laplacian matrix for the sinks) to solve this on a square domain, with Neumann boundary conditions at each side. The size of my domain is 100 x 100 points.

However, I wanted to know if there could be a faster way to solve for the steady state? I am solving for the field at each time step, which means that even if solver are relatively quick (a third of a second at the moment in c++ and Matlab), the time taken to solve for the field is a bottle neck for speeding up my code.

I have done some reading, and there appear to be some suggestions that Spectral methods could be used here, but I don't want to dive in only to find it takes longer. Similarly, some have suggested solving the top parabolic equation. Any advice here would be most appreciative.

Best,

Ben