clear all
f = @(y) [-2*y(2)*y(1);...
          y(1)^2+y(3)^2-y(2)^2-1;...
          -2*(y(2)+y(1))*y(3)];
jf = @(y) [-2*y(2),-2*y(1),0;...
           2*y(1),-2*y(2),2*y(3);...
           -2*y(3),-2*y(3),-2*(y(2)+y(1))];
tstar = 15;
y0 = [0.5;2;10];
ts = 1500;
k = tstar/ts;
tol = k/100;
y = zeros(3,ts+1);
y(:,1) = y0;
F = @(y,yn) y-k*f(y)-yn;
JF = @(y) eye(length(y))-k*jf(y);
for n = 1:ts
    y(:,n+1) = y(:,n);
    delta = -JF(y(:,n+1))\F(y(:,n+1),y(:,n));
    while (norm(delta) > tol)
        y(:,n+1) = y(:,n+1)+delta;
        delta = -JF(y(:,n+1))\F(y(:,n+1),y(:,n));
    end
end
t = linspace(0,tstar,ts+1);
E = (y(1,:).^2+y(2,:).^2+y(3,:).^2+1)./(2*y(1,:));
plot(t,E,'*')
