The equation for $u(x,t)$ is
$$
\partial_t u = \partial_{xx}u - u^3 + u + B\left(u_0-\frac{1}{L}\int_{0}^{L} u\;dx\right)
$$
This is a non-linear integro-differential equation (if not for the $u^3$, it would be linear). Laplace transform techniques don't get you very far with this equation, because the Laplace transform of $u^3$ is an integral in the Laplace domain. So even after transforming, you have an integro-differential equation; it doesn't simplify to an algebraic equation or ODE, as would happen with linear equations.
It won't be simple to develop your own, but numerical solutions are the way to go here. The authors of the paper from which this comes hastily mention that they use spatial discretization to reduce the problem to a system of ODEs in time, which they then solve with Runge-Kutta methods. The equation above is directly amenable to such an approach, but you have to formulate it. There are also a myriad of other techniques for numerical solutions to such equations. A good package I've been using recently is Chebfun for MATLAB, which can solve this using pseudospectral methods.
MATLAB and Chebfun
This is the easy way. The pde15s routine from the Chebfun package seems to solve this problem correctly. For a bump shaped initial condition with $L=4$, $B=10$, $u_0=1/2$, I get plots like

, which has the diffusion kind of properties I would expect from this equation. My numerical experiments tell me that in order to see interesting behavior, you have to have the equation on a ring with periodic boundary conditions (like the authors of that paper), or you need non constant initial conditions. For the constant initial conditions and Neumann boundary conditions like you propose, the solution doesn't even depend on $x$, i.e. the lack of initial $x$ variations can't be "smoothed out" any further by this kind of diffusion equation, so the solution is just a constant in $x$ along the rod that changes in time.
Finite Differences and Runge-Kutta
Another option is to do what the authors of that paper did. The idea is to sample $u(x,t)$ on a set of $N$ equidistant-spaced points (distance $\Delta x$), say $x_n\in[0,L]$, $n\in[1,2\ldots N]$. These are called nodes in a mesh. Let $u_n(t)=u(x_n,t)$ be the value of $u$ at the point $x_n$. You approximate your spatial derivatives and integrals with finite approximations, which in this case boils down to
$$
\partial_{xx} u(x_n,t)\approx \frac{u_{n+1}(t)-2u_{n}(t)+u_{n-1}(t)}{\left(\Delta x\right)^2}
$$
$$
\int_{0}^{L} u\;dx\approx \sum_{n=1}^N u_n(t) \Delta x
$$
The integral could be approximated with a higher order rule like the trapezoid or Simpson's rule. There is probably some optimum order for this problem (something like matching the approximation order used for the derivative), you'll have to find the details on how that works. Plugging these into your equation and enforcing the equation on the nodes gives a system of $N$ equations
$$
u_n' = \frac{u_{n+1}-2u_{n}+u_{n-1}}{\left(\Delta x\right)^2} - u_n^3 + u_n + B\left(u_0-\sum_{n=1}^N u_n \Delta x\right),\;\;\forall n\in[1,2\ldots N]
$$
I've dropped the explicit time dependence and let the prime denote a time derivative. These are $N$ ODEs in the $N$ unknowns, the $u_n(t)$ functions. You'll also need some more boundary conditions to specify what $u_{0}(t)$ and $u_{N+1}(t)$ are. Now you can use the Runge-Kutta method to solve this system of ODEs at discrete time steps. This is what the authors of the paper did.
You also have the option to directly discretize in space and time, which I think is the more standard approach. This is essentially a heat/diffusion type equation, so the discretization methods used there all apply. Check out the heat-equation example of finite differences on Wikipedia.