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