NDSolve 怎么提高精度

下面是代码:
n = 4;
Fig = Table[1, {i, 1, n}];
For[l = 1, l <= n, l++,
p1 = 1;
p2 = p1 (1 + 10^(-l - 13));
Clear[{s1, s2}];
s1 = y /.
NDSolve[{y'[x] == y[x] Cos[p1 x + y[x]], y[0] == 1},
y, {x, 0, 31}, AccuracyGoal -> 20, PrecisionGoal -> 20,
WorkingPrecision -> 64][[1]];
s2 = y /.
NDSolve[{y'[x] == y[x] Cos[p2 x + y[x]], y[0] == 1},
y, {x, 0, 31}, AccuracyGoal -> 20, PrecisionGoal -> 20,
WorkingPrecision -> 64][[1]];
Fig[[l]] = Plot[{s1[x] - s2[x]}, {x, 0, 30}, PlotRange -> Full];
]
Fig

试了很多方法, 但是, 只要p1 和p2 的区别小于10^(-17)左右, NDSolve 就算不出区别来,
于此比较, 例如p2 = p1 (1 + 10^-14) 就是比较清晰的, p2 = p1 (1 + 10^-15)就开始有细微震荡了.