clf, clear, hold on, format long, axis equal
% (x-x0)^2 + (y-y0)^2 = R^2

x = [1.1, 1.5, 2.2, 3.2, 4.5, 4.8]';
y = [0.7, 1.9, 3.1, 3.0, 2.3, 1.0]';

x0=0; y0=0; r0=2;
c = [x0, y0, r0]';

for iter=1:5
    r0m = ones(length(x),1) * r0;
    f = [(x - x0).^2 + (y - y0).^2 - r0^2];
    J = -2*[(x - x0), (y - y0), r0m];
    dc = -J\f;
    c = c + dc;
    x0 = c(1); y0 = c(2); r0 = c(3);
end
c
plot(x, y, '*');

phi = 0:pi/60:2*pi;
plot(c(1)+c(3)*cos(phi), c(2)+c(3)*sin(phi));