// Fonctions nécessaires pour le problème de Dirichlet


function y=f(x)
    y=(1+2*x-x.^2).*(exp(x));
endfunction

function y=c(x)
    y=x;
endfunction

// Résolution numérique du problème de Dirichlet


function U = dirichlet(f,c,a,b,h)

//Préparation de la discrétisation:

N = (1/h) - 1;
X = (0:h:1);

I = (2:N+1);

// Vecteur solution:
 
U = zeros(1,N+2);
U(1) = a;
U(N+2) = b;


// Matrice du problème discrétisé:

A = (1/(h*h))*(diag(2*ones(1,N)) + diag(-ones(1,N-1),-1) + diag(-ones(1,N-1),1)); //Laplacien standard

A = A + diag( c(X(I))) ; // Matrice "totale" avec le terme d'ordre 0:

// Résolution du problème matriciel

F = f(X(I)); //qui est un vecteur ligne!
F = F'; // Il devient une colonne!
F(1)= F(1) + (a/(h*h));
F(N) = F(N)+(b/(h*h));


V = A\(  F  )  // Résolution du système

U(I) = V' // Ma solution devient la partie "centrale" de U

// Tracons la solution 

//plot(X,U,'r-');
//plot(X,(1-X).*(exp(X)),'+');
//xtitle('Solution du problème de Dirichlet');
//

endfunction


function erreurlog(f,c,a,b)
    
// On prépare les graphiques d'erreur logarithmique
// On va comparer avec la fonction de référence (1-X)exp(X).
ERR = zeros(1,3);
H = 10.^(- (1:3));

// Ne pas se soucier de ca:
stacksize('max')


// On calcule l'erreur numérique (en norme infini) pour chaque valeur de h:
for i=1:3,
    h = H(i);
    U = dirichlet(f,c,a,b,h);
    X = (0:h:1);
    REF = (1-X).*(exp(X));
    ERR(i) = norm(U-REF,'inf');
end

//On trace maintenant les graphiques logarithmiques 

// Erreur numérique
plot2d('ll',H,ERR);

// Pente de référence
plot2d('ll',H,H.^2,style = 5);

endfunction
    
