clear all
A = 16*toeplitz([-2,1,0]);
y0 = [sqrt(2)/2;1;sqrt(2)/2];
f = @(y) A*y+y.^2;
tstar = 1;
tsrif = 1000;
ts = tsrif;
k = tstar/ts;
yERexp = y0;
for n = 1:ts
    Jn = A+diag(2*yERexp);
    yERexp = yERexp+k*phi1m(k*Jn)*f(yERexp);
end
yrif = yERexp;
tsrange = (10:10:100);
counter = 0;
for ts = tsrange
    counter = counter+1;
    k = tstar/ts;
    P = phi1m(k*A);
    yEexp = y0;
    yERexp = y0;
    for n = 1:ts
        yEexp = yEexp+k*P*f(yEexp);
        Jn = A+diag(2*yERexp);
        yERexp = yERexp+k*phi1m(k*Jn)*f(yERexp);
    end 
    errEexp(counter)=norm(yrif-yEexp,inf);
    errERexp(counter)=norm(yrif-yERexp,inf);
end
loglog(tsrange,errEexp,'*',tsrange,errERexp,'o')
  
