/// theta-schéma pour l'équation de diffusion u_t = d u_xx
/// en dimension 1 avec avec conditions au bord de Dirichlet 

clear();
stacksize(100000000);

d = 1; // coefficient de diffusion
T = 0.1; // temps final
L = 1; // taille du domaine [0;L]
M = 50; // nombre de mailles en espace
dx= L/(M+1); // pas en espace
x = (1:M)'*dx; // points de la grille 
xb = (0:M+1)'*dx; // grille incluant les bords
theta = 0.5; // choisir une valeur theta and [0;1]

if theta < 0.5 then
    dtmax = 1/(2*d*(1-theta))*dx*dx;
    rapport = 1/2;
    dt = min(rapport*dtmax,T);
    N = ceil(T/dt);
else
    N = 50;
    dt = T/N;
end

nb = d*dt/(dx^2);
A = sparse(-2*eye(M,M)+diag(ones(M-1,1),1)+diag(ones(M-1,1),-1));
CI = sin(2.*%pi*x);
CIb = sin(2.*%pi*xb);
B = sparse(eye(M,M)-theta*nb*A);
R = sparse(eye(M,M)+(1-theta)*nb*A);
[B,rk] = lufact(B);

u = CI;
t = 0;
sol = zeros(M+2,N+1);
sol(:,1) = CIb;

for n=0:N-1
    t = t+dt;
    rhs = R*u;
    u = lusolve(B,rhs);
    ub = [0;u;0];
    sol(:,n+2) = ub;
    drawlater();
    clf();
    a=gca(); 
    a.data_bounds=[0,-1;1,1];
    plot(xb,ub,'d',xb,exp(-4.*%pi^2*t)*CIb);
    xlabel("x")
    drawnow();
end

halt();

clf();
subplot(211)
plot(xb,[0;u;0],'d',xb,exp(-4.*%pi^2*t)*CIb);
subplot(212)
plot(x,abs(u-exp(-4.*%pi^2*t)*CI))

halt();
clf();
surf(0:dt:T,xb,sol)
xlabel("t")
ylabel("x")
