下面是代码:
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)就开始有细微震荡了.

嗬，我好像还是第一次在吧里见到会调WorkingPrecision并且还知道针对性地调整AccuracyGoal与PrecisionGoal的人，虽然说要是能不用For就更完美了吧。
答案很简单：Plot它也有WorkingPrecision：
n = 4;
Fig = Table[p1 = 1;
p2 = p1 (1 + 10^(-l - 13));
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]];
Plot[{s1[x] - s2[x]}, {x, 0, 30}, PlotRange -> All, WorkingPrecision -> 64], {l,
n}]