clear all
close all
c = 0.8;
d = 0.01;
rho = 50;
tstar = 0.1;
m = 100;
h = 1 / (m-1);
Pe = abs(c)*h/(2*d)
x = linspace(0,1,m)';
D1 = toeplitz(sparse(1,2,-1/(2*h),1,m),sparse(1,2,1/(2*h),1,m));
D1(m,m-1:m) = [0,0];
D2 = toeplitz(sparse([1,1],[1,2],[-2,1]/h^2,1,m));
D2(m,m-1:m) = [2,-2]/h^2;
A = d*D2-c*D1;
A(1,[1:2]) = [0,0];
g = @(y) rho*y.*(y-1/2).*(1-y);
dg = @(y) rho*((y-1/2).*(1-y)+y.*(1-y)-y.*(y-1/2));
b = @(t,y) [0;g(y(2:m))];
jb = @(t,y) diag([0;dg(y(2:m))]);
y0 = 5*x.*(1-x).^2;
tsrif = 1000;
ts = tsrif;
k = tstar/ts;
y = y0;
t = 0;
tol = min(h^2,k)/100;
for n=1:ts
  yn = y;
  F = @(y) y-yn-k*(A*y+b(t,y));
  JF = @(y) speye(length(y))-k*A-k*jb(t,y);
  y = yn;
  delta = -JF(y)\F(y);
  while (norm(delta,inf) > tol)
    y = y+delta;
    delta = -JF(y)\F(y);
  end
  y = y+delta;
  t = t+k;
end
yrif = y;
tsrange = (10:10:50);
counter = 0;
for ts = tsrange
  counter = counter + 1;
  k = tstar/ts;
  y = y0;
  t = 0;
  tol = min(h^2,k)/100;
  for n=1:ts
    yn = y;
    F = @(y) y-yn-k*(A*y+b(t,y));
    JF = @(y) speye(length(y))-k*A-k*jb(t,y);
    y = yn;
    delta = -JF(y)\F(y);
    while (norm(delta,inf) > tol)
      y = y+delta;
      delta = -JF(y)\F(y);
    end
    y = y+delta;
    t = t+k;
  end 
  err(counter) = norm(y-yrif,inf);
end
loglog(tsrange,err,'*',tsrange,tsrange.^(-1))
