// discrétisation en temps/espace
N = 100;
a = 0
b = 1;
dx = (b-a)/(N-1);
x = a:dx:b;

nu = 1;
dt = 0.1;
// matrice A du Laplacien avec conditions périodiques au bord
e = ones(N,1);
ee = ones (N-1,1);
A = diag(ee,-1) + diag(ee,+1) - 2*diag(e);

// Change boundary conditions : 
//-------------------------------
// a) Dirichlet boundary conditions

// b) periodic boundary conditions
//A(1,N) = 1;
//A(N,1) = 1;
// c) Neumann boundary conditions
//A(1,1) = -1
//A(N,N) = -1 
// rescale with respect to lambda=dt/dx^2 
A = A*(dt/dx^2);

B = sparse(eye(N,N)-nu*A);

u0 = 0.5*(exp(-(x+4).^2)+exp(-(x-4).^2))';
u = u0;
t = 0;
Tf = 50;

while t < Tf
    t = t+dt;
    //u1 = B\u;
    u1 = ( eye(N,N) + A )*u 
    drawlater();
    clf();
    plot(x,u','r')
    drawnow();
    u = u1;
end





