Instability of FindRoot solution

Using FindRoot on a smooth function of a parameter and a variable, the resulting solution should be a smooth function of the parameter, by the Implicit Function Theorem. However, Mathematica returns a very spiky unstable function. One way to fix the problem seems to avoid Newton’s method. Are there others, in particular that work with more than one variable?

MWE:

Clear[cdf, s0, s, sol] cdf[x_] := -2 (x - Sinh[1/2] +       Sinh[1/2 - x])/(-2 (1 - Sinh[1/2] + Sinh[1/2 - 1])) s0[a_] :=   s0[a] = x /. FindRoot[D[x*(cdf[a - x] - cdf[x]), x] == 0, {x, 0.5}] Plot[s0[a], {a, 0, 1}] (*A different root finding method does not cause instability.*) s[a_] := s[a] =    x /. FindRoot[D[x*(cdf[a - x] - cdf[x]), x] == 0, {x, 0.01, 0.99}] Plot[s[a], {a, 0, 1}] (*The problem is worse with two variables*) sol[c_] :=   sol[c] = {p1, p2} /.     FindRoot[{D[p1*(cdf[p2 - p1] - cdf[p1]), p1] == 0,   D[(p2 - c)*(1 - cdf[p2 - p1]), p2] == 0}, {{p1, 0.2}, {p2, 0.5}}] Plot[{sol[c][[1]], sol[c][[2]]}, {c, 0, 1}] (*Using non-Newton methods of root finding does not help with two variables.*) sol[c_] :=   sol[c] = {p1, p2} /.     FindRoot[{D[p1*(cdf[p2 - p1] - cdf[p1]), p1] == 0,       D[(p2 - c)*(1 - cdf[p2 - p1]), p2] == 0}, {{p1, 0.1, 0.8}, {p2,        0.2, 0.9}}] Plot[{sol[c][[1]], sol[c][[2]]}, {c, 0, 1}]