function [err, grad] = lserr(p, x, y, fun)

err = sum(y - fun(p(1), p(2), x)).^2;

if nargout > 1
    grad = lserr_gradient(p, x, y);
end