clear all
close all
m = 3;
A = (m+1)^2*toeplitz([-2,1,zeros(1,m-2)]);
b = @(y) y.^2;
f = @(y) A*y+b(y);
jf = @(y) A+2*diag(y);
tstar = 1;
ts = 10;
k = tstar/ts;  
yee = zeros(m,ts+1);
yei = zeros(m,ts+1);
yeexp = zeros(m,ts+1);
tol = k/100;
F = @(y,yn) y-yn-k*f(y);
JF = @(y) eye(m)-k*jf(y);
y0 = sin((1:m)*pi/(m+1));
yee(:,1) = y0;
yei(:,1) = y0;
yeexp(:,1) = y0;
P = phi1m(k*A);
for n = 1:ts
  yee(:,n+1) = (eye(m)+k*A)*yee(:,n)+k*b(yee(:,n));
  yei(:,n+1) = yei(:,n);
  delta = -JF(yei(:,n+1))\F(yei(:,n+1),yei(:,n));
  while (norm(delta,inf) > tol)
    yei(:,n+1) = yei(:,n+1)+delta;
    delta = -JF(yei(:,n+1))\F(yei(:,n+1),yei(:,n));
  end
  yeexp(:,n+1) = yeexp(:,n)+k*P*f(yeexp(:,n));
end
t = linspace(0,tstar,ts+1);
plot(t,yee)
figure
plot(t,yei)
figure
plot(t,yeexp)