clear all
close all
m = 100;
x = linspace(0,1,m)';
h = 1/(m-1);
A = toeplitz(sparse([1,1],[1,2],[-2,1]/h^2,1,m));
B = toeplitz(sparse(1,2,-1/(2*h),1,m),sparse(1,2,1/(2*h),1,m));
A(m,m-1:m) = [2,-2]/h^2; % Neumann
B(m,m-1:m) = [0,0]; % Neumann
d = 0.01;
c = 0.8;
rho = 50;
C = d*A-c*B;
C(1,1:2) = [0,0]; % Dirichlet
%C(m,m-1:m) = [0,0]; % Diriclet
Peclet = c*h/(2*d)
g = @(u) rho*u.*(u-1/2).*(1-u);
tstar = 1;
timesteps = 100;
k = tstar/timesteps;
y0 = 5*x.*(1-x).^2;
plot(x,y0)
axis([0,1,-0.1,1])
pause
y = y0;
P = phi1m(k*C);
for n = 1:timesteps
  y = y+k*P*(C*y+[0;g(y(2:m))]);
  plot(x,y)
  axis([0,1,-0.1,1])
  pause(0.1)
end