Solving a PDE gives stiffness error

Question

I'm solving a PDE with

Needs["DifferentialEquations`NDSolveProblems`"];
Needs["DifferentialEquations`NDSolveUtilities`"];

tmin = -2; tmax = 2;
(* u[x, 0] \[Equal] 4Sech[x]^2*)

ν = 1; ϵ = 0.001;(*0.005*);

xmin = -2; xmax = 2; 
sol = 
 NDSolve[ {D[u[x, t], t] + 
     u[x, t]*D[u[x, t], x] + ϵ^2*D[u[x, t], {x, 3}] == ν*
     D[u[x, t], {x, 2}], u[x, 0] == Cos[x], u[xmin, t] ==  u[xmax, t] }, 
  u, {x, xmin, xmax}, {t, tmin, tmax}, Method -> "Automatic"]

During evaluation of In[40]:= NDSolve::mxsst: Using maximum number of grid points 10000 allowed by the MaxPoints or MinStepSize options for independent variable x.

During evaluation of In[40]:= NDSolve::ndsz: At t == -5.72262*10^-7, step size is effectively zero; singularity or stiff system suspected.

During evaluation of In[40]:= NDSolve::eerr: Warning: scaled local spatial error estimate of 5.058814843192396`^10 at t = -5.7226210^-7 in the direction of independent variable x is much greater than the prescribed error tolerance. Grid spacing with 10001 points may be too large to achieve the desired accuracy or precision. A singularity may have formed or a smaller grid spacing can be specified using the MaxStepSize or MinPoints method options.

(* Out[40]= {{u -> -7 InterpolatingFunction[{{…, -2., 2., …}, {-5.72262 10 , 2.}}, <>]}} *)

Solve Burgers in form ν ∂u^2/∂^2x= -(∂u/∂t)- u ∂u/∂x-ϵ^2 ∂^3u/∂x^3

Plot3D[u[x, t] /. Flatten[sol], {x, xmin, xmax}, {t, tmin, tmax}, 
 PlotPoints -> 100, PlotRange -> All, ColorFunction -> "Rainbow"]

But can't get an output because

NDSolve::ndsz: At t == -5.72262*10^-7, step size is effectively zero; singularity or stiff system suspected.

Any ideas how to fix?

Answer 1

Note: in the comment section under the OP we found a typo, which has been taken into consideration in the suggested solution below.

The following runs without any errors and complaints

tmin = -2;
tmax = 2;
ν = 1;
ε = 0.001;
xmin = -2;
xmax = 2;
points = 2000;

sol = NDSolveValue[{D[u[x, t], t] + 
     u[x, t]*D[u[x, t], x] + ε^2*D[u[x, t], {x, 3}] == ν*
     D[u[x, t], {x, 2}], u[x, -2] == Cos[x], 
   u[xmin, t] == u[xmax, t]}, 
  u[x, t], {x, xmin, xmax}, {t, tmin, tmax}, 
  Method -> {"MethodOfLines", 
    "SpatialDiscretization" -> {"TensorProductGrid", 
      "MaxPoints" -> points, "MinPoints" -> points, 
      "DifferenceOrder" -> "Pseudospectral"}, Method -> "LSODA"}, 
  MaxSteps -> Infinity]

Visualizing

Plot3D[sol, {x, xmin, xmax}, {t, tmin, tmax}]