Numerical instabilities in a solution to a partial differential equation

I have been trying to solve a partial differential equation known as KdV (Korteweg-de Vries) equation using NDSolve.
KdV equation is written as:

$frac{partial u(x,t)}{partial t}+frac{partial^3 u(x,t)}{partial x^3}=-6u(x,t)frac{partial u(x,t)}{partial x}$

One of the analytical solution of this equation is a Jacobi elliptical function
$u_{A} = text{A JacobiCN}[B*(x- c t),m]^2$
where A, B, c and m are known parameters.

So the initial condition is
$u_0 = text{A JacobiCN}[B*x,m]^2 $.

The boundary condition is: $u(t,-a) = u(t,a)$.

Code for analytical solution is

a = 30; m = 0.9; A = 0.9;
B = Sqrt[1.5*A] / Sqrt[2*m + m^2]
c = 6*(-A + A*m + A*m^2)/(m*(2 + m))
uA[x_,t_] := A*JacobiCN[B*(x - c*t), m]^2;
Plot[uA[x,0.5], {x, -a, a}, LabelStyle ->  Directive[14, Bold,   Black], AxesLabel -> {Style["x", 14, Bold, Black], Style["uA", 14, Bold, Black]}] 

Plot of analytical solution for t = 0.5

When I solve the partial differential equation using NDSolve, the numerical instabilities arise and the numerical solution doesn’t match the analytical solution.

eq = {D[u[x, t], {t, 1}] + D[u[x, t], {x, 3}] == -6 u[x, t] D[u[x, t], {x, 1}], u[x, 0] == A*JacobiCN[B*x, m]^2, u[-a, t] == u[a, t]};
sol = Flatten@NDSolve[eq, u, {t, 0, 1}, {x, -a, a}]
Plot[Evaluate[u[x, 0.5] /. sol], {x, -a, a}, LabelStyle ->  Directive[14, Bold, Black], AxesLabel -> {Style["x", 14, Bold, Black], Style["u", 14, Bold, Black]}]

I am also getting the warning message
NDSolve::mxsst: Using maximum number of grid points 10000 allowed by the MaxPoints or MinStepSize options for independent variable x.
I have tried different options and methods given in NDSolve documentation and also tried solutions given in https://mathematica.stackexchange.com/questions/9277/i-failed-to-solve-a-set-of-one-dimension-fluid-mechanics-pdes-with-ndsolve/9291#9291 but nothing seems to be working.
Following method removed the warning message but there was no improvement in the numerical solution.

nxy = 500;
sol = Flatten@NDSolve[eq, u, {t, 0, 1}, {x, -a, a}, Method -> {"MethodOfLines", 
     "SpatialDiscretization" -> {"TensorProductGrid", "MaxPoints" -> nxy, "MinPoints" -> nxy, 
       "DifferenceOrder" -> "Pseudospectral"}, Method -> "Adams"}, MaxSteps -> [Infinity]]