clear all
close all
tstar = 1;
f = @(t,y) y*(1-y);
df = @(t,y) (1-y)-y;
y0 = 1/2;
theta = 1/2;
tsrif = 1000;
ts = tsrif;
k = tstar/ts;
tol = k^2/100;
F = @(tn,y,yn) y-yn-k*(1-theta)*f(tn,yn)-k*theta*f(tn+k,y); 
J = @(tn,y) eye(length(y))-k*theta*df(tn+k,y);
t = linspace(0,tstar,ts+1);
y = zeros(length(y0),ts+1);
y(:,1) = y0;
for n = 1:ts
    y(:,n+1) = y(:,n);
    delta = -J(t(n),y(:,n+1))\F(t(n),y(:,n+1),y(:,n));
    while (norm(delta,inf) > tol)
        y(:,n+1) = y(:,n+1)+delta;
        delta = -J(t(n),y(:,n+1))\F(t(n),y(:,n+1),y(:,n));
    end
end
yrif = y(:,ts+1);
tsrange = (10:10:100);
counter = 0;
for ts = tsrange
    counter = counter+1;
    k = tstar/ts;
    yn = y0;
    tol = k^2/100;
    F = @(tn,y,yn) y-yn-k*(1-theta)*f(tn,yn)-k*theta*f(tn+k,y); 
    J = @(tn,y) eye(length(y))-k*theta*df(tn+k,y);  
    t = linspace(0,tstar,ts+1);
    for n = 1:ts
        y = yn;
        delta = -J(t(n),y)\F(t(n),y,yn);
        while (norm(delta,inf) > tol)
            y = y+delta;
            delta = -J(t(n),y)\F(t(n),y,yn);  
        end
        yn = y;
    end
    err(counter) = norm(y-yrif,inf);
end
loglog(tsrange,err,'*',tsrange,err(end)*(tsrange/tsrange(end)).^(-2))
