clear all
close 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)];
df = @(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))];
y0 = [0.5;2;10];
tstar = 15;
tsrif = 5000;
ts = tsrif;
%
k = tstar/ts;
tol = k^2/100;
y = zeros(3,ts+1);
y(:,1) = y0;
F = @(yn,x) x-yn-k*f(x);
JF = @(x) eye(3)-k*df(x);
n = 1;
y(:,n+1)=y(:,n);
delta = -JF(y(:,n+1))\F(y(:,n),y(:,n+1));
while (norm(delta,inf)>tol)
    y(:,n+1) = y(:,n+1)+delta;
    delta = -JF(y(:,n+1))\F(y(:,n),y(:,n+1));
end
y(:,n+1) = y(:,n+1)+delta;
F = @(yn,yn1,x) x-4*yn1/3+yn/3-k*2*f(x)/3;
JF = @(x) eye(3)-k*2*df(x)/3;
for n = 1:ts-1
    y(:,n+2)=y(:,n+1);
    delta = -JF(y(:,n+2))\F(y(:,n),y(:,n+1),y(:,n+2));
    while (norm(delta,inf)>tol)
        y(:,n+2) = y(:,n+2)+delta;
        delta = -JF(y(:,n+2))\F(y(:,n),y(:,n+1),y(:,n+2));
    end
    y(:,n+2) = y(:,n+2)+delta;
end
yrif = y;
E = (y(1,:).^2+y(2,:).^2+y(3,:).^2+1)./(2*y(1,:));
t = linspace(0,tstar,ts+1);
plot(t,E,'*')
tsrange = 1000:100:1500;
counter = 0;
for ts = tsrange
    counter = counter+1;
    k = tstar/ts;
    tol = k^2/100;
    y = zeros(3,ts+1);
    y(:,1) = y0;
    F = @(yn,x) x-yn-k*f(x);
    JF = @(x) eye(3)-k*df(x);
    n = 1;
    y(:,n+1)=y(:,n);
    delta = -JF(y(:,n+1))\F(y(:,n),y(:,n+1));
    while (norm(delta,inf)>tol)
        y(:,n+1) = y(:,n+1)+delta;
        delta = -JF(y(:,n+1))\F(y(:,n),y(:,n+1));
    end
    y(:,n+1) = y(:,n+1)+delta;
    F = @(yn,yn1,x) x-4*yn1/3+yn/3-k*2*f(x)/3;
    JF = @(x) eye(3)-k*2*df(x)/3;
    for n = 1:ts-1
        y(:,n+2)=y(:,n+1);
        delta = -JF(y(:,n+2))\F(y(:,n),y(:,n+1),y(:,n+2));
        while (norm(delta,inf)>tol)
            y(:,n+2) = y(:,n+2)+delta;
            delta = -JF(y(:,n+2))\F(y(:,n),y(:,n+1),y(:,n+2));
        end
        y(:,n+2) = y(:,n+2)+delta;
    end        
    err(counter) = norm(y(:,ts+1)-yrif(:,tsrif+1),inf);
end
figure
loglog(tsrange,err,'*',...
       tsrange,err(end)*(tsrange/tsrange(end)).^-2)
