function K = krivulja(f, gradf, T0, delta, n)
%K = krivulja(f, gradf, T0, delta, n) returns a sequence of points on
%the implicitly given curve f(x, y) = 0.
%f ....... the function (of two variables) f,
%gradf ... the gradient of f,
%T0 ...... starting point on the curve,
%delta ... Euclidean distance between two consecutive points,
%n ....... number of points.
%system II (correction to T0 - to obtain a point on the curve)
g = gradf(T0(1), T0(2));
ng = [-g(2); g(1)];
%define the function (and its Jacobi matrix) corresponding to system II
F = @(X) [f(X(1), X(2)); ng'*(X - T0)];
JF = @(X) [gradf(X(1), X(2))'; ng'];
%use T0 as an initial guess and solve system II using Newton's method
K(:, 1) = newton(F, JF, T0);
%system I (obtaining succesive points on the curve)
for k = 2:n
%define the function (and its Jacobi matrix) corresponding to system I
F = @(X) [f(X(1), X(2)); (X - K(:, k-1))'*(X - K(:, k-1)) - delta^2];
JF = @(X) [gradf(X(1), X(2))'; 2*(X - K(:, k-1))'];
%determine the tangent through current point
g = gradf(K(1,k-1), K(2,k-1));
ng = [-g(2); g(1)];
%evaluate the next point
K(:, k) = newton(F, JF, K(:, k-1) + delta*ng/norm(ng));
%break the loop once close enough to the starting point (in case of closed curve)
if(norm(K(:, k) - K(:, 1)) < delta-10*eps)
break;
end
end