clear all
tstar = 1;
mrange = [10:10:100];
counter = 0;
for m = mrange
    counter = counter+1;
    k = tstar/m;
    y = zeros(1,m+1);
    y0 = 1/2;
    y(1) = y0;
    f = @(t,y) y.*(1-y);
    jf = @(t,y) 1-y-y;
    F = @(tn,yn,t,y) y-yn-k/2*(f(tn,yn)+f(t,y));
    JF = @(t,y) 1-k/2*jf(t,y);
    t = 0;
    tol = 1e-12;
    for n = 1:m
        y(n+1) = y(n);
        delta = -JF(t+k,y(n+1))\F(t,y(n),t+k,y(n+1));
        while (norm(delta,inf) > tol)
            y(n+1) = y(n+1)+delta;
            delta = -JF(t+k,y(n+1))\F(t,y(n),t+k,y(n+1));
        end
        t = t+k;
    end
    t = linspace(0,tstar,m+1);
    errore(counter) = norm(y-exp(t)./(1+exp(t)),inf);
end
loglog(mrange,errore,'*',mrange,errore(1)*(mrange/mrange(1)).^(-2))