UpdateRemark
Even further check shows
NDSolvehas used 2-norm instead of infinity-norm for priori error estimates since v9, see this post for more info.
Increasing points to 2000 doesn't significantly change the solution, soand the result is consistent with that of v5.2:
So I guess the solution is reliable.
I tested the sample in other versions, and found the backslide is introduced in v9:
v9.0.1 
v8.0.4 
Further check shows NDSolve has used 9615 grid points in v5.2, and 9681 in v8.0.4 for spatial discretization, thus this seems to be a backslide of priori error estimates. If we set the Method -> {"MethodOfLines", "SpatialDiscretization" -> {"TensorProductGrid", "MaxPoints" -> 9615, "MinPoints" -> 9615, "DifferenceOrder" -> 4}, Method -> "StiffnessSwitching"} in higher versions:
NDSolve gives the desired result. Notice Method -> "StiffnessSwitching" is necessary here.
But still, my method that takes the consevation law into consideration is cheaper. (Only 1500 grid points are needed. )
Increasing points to 2000 doesn't significantly change the solution, so I guess the solution is reliable.
Increasing points to 2000 doesn't significantly change the solution, and the result is consistent with that of v5.2:
So I guess the solution is reliable.
I tested the sample in other versions, and found the backslide is introduced in v9:
v9.0.1 v8.0.4
Further check shows NDSolve has used 9615 grid points in v5.2, and 9681 in v8.0.4 for spatial discretization, thus this seems to be a backslide of priori error estimates. If we set the Method -> {"MethodOfLines", "SpatialDiscretization" -> {"TensorProductGrid", "MaxPoints" -> 9615, "MinPoints" -> 9615, "DifferenceOrder" -> 4}, Method -> "StiffnessSwitching"} in higher versions:
NDSolve gives the desired result. Notice Method -> "StiffnessSwitching" is necessary here.
But still, my method that takes the consevation law into consideration is cheaper. (Only 1500 grid points are needed. )