clear all
close all
% esercizio 4 PDE, differenze finite, Eulero implicito, tstar = 0.1
a = 0;
b = 1;
m = 50;
x = linspace(a,b,m)';
u0 = 5*x.*(1-x).^2;
h = (b-a)/(m-1);
A = toeplitz(sparse([1,1],[1,2],[-2,1]/h^2,1,m));
A(m,m-1:m) = [2,-2]/h^2; % Neumann
%A(m,m-1:m) = [0,0]; % Dirichlet
B = toeplitz(sparse(1,2,-1/(2*h),1,m),...
    sparse(1,2,1/(2*h),1,m));
B(m,m-1:m) = [0,0]; % Neumann, Dirichlet
c = 8;
d = 0.1;
Peclet = c*h/(2*d)
rho = 10;
C = d*A-c*B;
C(1,1:2) = [0,0];
g = @(u) rho*u.*(u-1/2).*(1-u);
dg = @(u) rho*((u-1/2).*(1-u)+u.*(1-u)-u.*(u-1/2));
tstar = 0.1;
tsrif = 1000;
ts = tsrif;
k = tstar/ts;
F = @(y,yn) y-yn-k*C*y-k*[0;g(y(2:m))]; % Neumann
%F = @(y,yn) y-yn-k*C*y-k*[0;g(y(2:m-1));0]; % Dirichlet
J = @(y) speye(length(y))-k*C-k*spdiags([0;dg(y(2:m))],0,m,m); % Neumann
%J = @(y) speye(length(y))-k*C-k*spdiags([0;dg(y(2:m-1));0],0,m,m); % Dirichlet
yn = u0;
tol = 1e-6;
for n = 1:ts
    y = yn;
    delta = -J(y)\F(y,yn);
    while (norm(delta,inf) > tol)
        y = y+delta;
        delta = -J(y)\F(y,yn);
    end
    plot(x,y)
    axis([0,1,-0.2,1.2])
    pause(0.1)
    yn = y;
end
yrif = y;
tsrange = 100:10:200;
counter = 0;
for ts = tsrange
  counter = counter+1;
  k = tstar/ts;
  F = @(y,yn) y-yn-k*C*y-k*[0;g(y(2:m))];
  J = @(y) speye(length(y))-k*C-k*spdiags([0;dg(y(2:m))],0,m,m);
  yn = u0;
  tol = 1e-6;
  for n = 1:ts
    y = yn;
    delta = -J(y)\F(y,yn);
    while (norm(delta,inf) > tol)
      y = y+delta;
      delta = -J(y)\F(y,yn);
    end
    yn = y;
  end
  err(counter) = norm(y-yrif,inf);
end
loglog(tsrange,err,'*',tsrange,err(end)*(tsrange/tsrange(end)).^(-1))
