Another problem related to conservation law. Based on the experience obtained in e.g. here, let's avoid symbolic expansion of D. I'll use pdetoode for the task.
q[x_] = (Erf[x] - 1)/2 - 5 Sech[x - 1]; xmax = 300; tmax = 15; With[{u = u[x, t], mid = mid[x, t]}, eq = {D[u, t] == D[mid, x], mid == -D[u, x, x] + 3 u^2}]; ic = u[x, 0] == q[x]; points = 1500; difforder = 2; domain = {-xmax, xmax}; grid = Array[# &, points, domain]; (* Definition of pdetoode isn't included in this post, please find it in the link above. *) ptoofunc = pdetoode[{u, mid}[x, t], t, grid, difforder]; odeadd = ptoofunc /@ eq[[2]]; ode = Block[{mid}, Set @@ odeadd; eq[[1]] // ptoofunc]; odeic = ptoofunc@ic; sollst = NDSolveValue[{ode, odeic}, u /@ grid, {t, 0, tmax}]; // AbsoluteTiming sol = rebuild[sollst, grid, 2]; // AbsoluteTiming rstlst = Plot[sol[x, #] // Evaluate, {x, -xmax, xmax}, PlotRange -> {-8, 1}, PlotPoints -> 50] & /@ Range[0, tmax, 0.5]; // AbsoluteTiming ListAnimate[rstlst, ControlPlacement -> Top] 
Oh, I didn't set the boundary conditions, but given the solution is localized i.e. there's no interaction between boundary and solution, this should not be too big a problem.
Increasing points to 2000 doesn't significantly change the solution, so I guess the solution is reliable.