The easiest solution really depends on the initial conditions and the coefficients. If your coefficients are not too large, quite frankly the easiest way to do so would be a Method of Lines. Approximate the LaPlacian to order 2 in space by the tridiagonal matrix (1,-2,1) and now you have a system of ODEs. In MATLAB if you define A like that then if you loop
u(t+1) = alpha/k * u(t) .*(1-u(t).^q) + A*u(t);
with a small enough $\Delta t$ you will get a solution. If you're just looking for answers, I'd try this first. If that's close but you need it slightly better, plug it into ode45. If it's stiffer than that, ode15s.
That shouldn't take more than an hour and if none of those are sufficient, then I would suggest doing something fancy. As Magneto suggested a fully implicit Crank-Nicholson will get the job done since it's order(2,2) and unconditionally stable, but you'll need to solve implicit equations (via fsolve). Spectral solutions that Mabuza suggested also should work well. If are solving a very stiff version and have a GPGPU using MATLAB's GPU FFTs for a spectral solution will probably be the best answer (in terms of "computational power per man hours").
But seriously, always try the easy way first!