syms a b x y r
eqn = (x-a)^2 + (y-b)^2 - r^2
J = jacobian(eqn, [x y])
syms tolp toln
assume(tolp > 0 & toln > 0)
range_error = 0.5
eqn1 = eqn + tolp - range_error^2
eqn2 = eqn - toln - range_error^2
J1 = jacobian(eqn1, [x y])
J2 = jacobian(eqn2, [x y])
What I did was transform
%sqrt((x-a)^2 + (y-b)^2) = r + delta, -range_error <= delta <= range_error
into equalities, giving a name to the difference between the ideal match and the actual match; split it into two parts, one with a positive difference and one with a negative difference, and require that the variable be positive. Like A > B means that A = B + delta where delta > 0. MATLAB is a lot more comfortable reasoning about equalities and then eliminating the branches that would violate the assume(), than it is trying to solve inequalities.
