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