function y = at(f, xk)
    syms x y;
    y = double(subs(f, [x, y], [xk(1), xk(2)]));
end