clear all
close all
theta = 1/2; % trapezi
if (theta == 1/2)
    order = 2;
else
    order = 1;
end
mrif = 500;
m = mrif;
tstar = 1;
y0 = 1/2;
f = @(t,y) y.*(1-y);
jf = @(t,y) (1-y)-y;
k = tstar/m;
F = @(t,y,tn,yn) y-k*theta*f(t,y)-...
    yn-k*(1-theta)*f(tn,yn);
J = @(t,y) eye(length(y))-k*theta*jf(t,y);
y = zeros(length(y0),m+1);
y(:,1) = y0;
tn = 0;
tol = k^2/100;
for n = 1:m
    y(:,n+1) = y(:,n);
    delta = -J(tn+k,y(:,n+1))\F(tn+k,y(:,n+1),tn,y(:,n));
    while (norm(delta,inf) > tol)
        y(:,n+1) = y(:,n+1)+delta;
        delta = -J(tn+k,y(:,n+1))\F(tn+k,y(:,n+1),tn,y(:,n));
    end
    tn = tn+k;
end
yrif = y;
plot(linspace(0,tstar,m+1),yrif,'*')
mrange = (10:20:100);
counter = 0;
for m = mrange
    counter = counter+1;
    k = tstar/m;
    F = @(t,y,tn,yn) y-k*theta*f(t,y)-...
        yn-k*(1-theta)*f(tn,yn);
    J = @(t,y) eye(length(y))-k*theta*jf(t,y);
    y = zeros(length(y0),m+1);
    y(:,1) = y0;
    tn = 0;
    tol = k^2/100;
    for n = 1:m
        y(:,n+1) = y(:,n);
        delta = -J(tn+k,y(:,n+1))\F(tn+k,y(:,n+1),tn,y(:,n));
        while (norm(delta,inf) > tol)
            y(:,n+1) = y(:,n+1)+delta;
            delta = -J(tn+k,y(:,n+1))\F(tn+k,y(:,n+1),tn,y(:,n));
        end
        tn = tn+k;
    end
    err(counter) = norm(y(:,m+1)-yrif(:,mrif+1),inf);
end
figure
loglog(mrange,err,'*',mrange,err(1)*(mrange/mrange(1)).^(-order))