clear all
close all
m = 101;
a = 0;
b = 1;
tstar = 1;
x = linspace(a,b,m)';
h = (b-a)/(m-1);
y0 = 10*x.^2.*(1-x).^2;
A = toeplitz(sparse([1,1],[1,2],[-2/h^2,1/h^2],1,m)/100);
% Neumann conditions
A(1,1:2) = [-2/h^2,2/h^2]/100;
A(m,m-1:m) = [2/h^2,-2/h^2]/100;
g = @(u) 10*u.*(u-1/2).*(1-u);
dg = @(u) 10*((u-1/2).*(1-u)+u.*(1-u)-u.*(u-1/2));
a11 = 0;
a21 = 1/2;
a22 = 1/2;
b1 = 1/2;
b2 = 1/2;
tol = 1e-5;
ts = 100;
k = tstar/ts;
y = y0;
for n = 1:ts
% f1
  f1 = A*y+g(y);
% f2  
  F = @(f2) f2-A*y-A*k*a21*f1-A*k*a22*f2-g(y+k*a21*f1+k*a22*f2);
  J = @(f2) speye(size(A))-A*k*a22-k*a22*diag(dg(y+k*a21*f1+k*a22*f2));
  f2 = A*y+g(y);
  errest = -J(f2)\F(f2);
  while (norm(errest,inf) > tol)
    f2 = f2+errest;
    errest = -J(f2)\F(f2);
  end
% y_{n+1}
  y = y+k*b1*f1+k*b2*f2;
  plot(x,y)
  axis([a,b,0,1])
  pause(0.1)
end
