////---------------------------------------------------
////
//// Resolution du Laplacien 1 D en differences finies.
////
////---------------------------------------------------
//// -u"=f sur [a,b]
//// condition de Dirichlet en a: u(a)=alpha u(b)=beta.
//// condition de Neumann en b: u'(b)=g
//// Auteur : V. Martin
//=====================================================
choix=2
//===========================
// Parametres physiques
a = 0;                           // taille du domaine
b = %pi;
alpha = 0;                       // conditions de Dirichlet
deff('y=f(x)','y=-((12*(x-1).^2).*exp(x)+8*((x-1).^3).*exp(x)+((x-1).^4).*exp(x))')
deff('y=solex(x)','y=((x-1).^4).*exp(x)')
x0=0;
g=4*(x0-1).^3.*exp(x0)+(x0-1).^4.*exp(x0)
beta = solex(b)

//===========================
// Parametres numeriques
n = 10                       // nb de points de discretisation
E=[]
H=[]
for i=1:10

    h = (b-a)/(n+1)
    
    //===========================
    // Maillage
    x = linspace(a,b,n+2)
    
    //===========================
    // Assemblage de la matrice
    e0 = ones(n+1,1)
    e1 = ones(n,1)
    A = 2*diag(sparse(e0),0)-diag(sparse(e1),-1)-diag(sparse(e1),1)
    A = A/h/h
    
    if choix==1 then
        A(1,1)=-1/h
        A(1,2)=1/h
    end
    if choix==2 then
        A(1,1:3)=[-3/2/h,2/h,-1/2/h]
    end
    
    
    //===========================
    // Second membre
    
    F = f(x(1:$-1))
    F(1) = g
    F(n+1) = F(n+1)+beta/h/h
    
    //===========================
    // Resolution
    u = A\F'
    u=[u;beta]
    
    //===========================
    // Calcul de l'erreur
    E=[E;max(abs(u-solex(x')))]
    H=[H;h]

    n=n*2    
end
plot(log10(H),log10(E))
plot(log10(H),1+log10(H),'r-.')
plot(log10(H),1+2*log10(H),'m*-')
legend('erreur du schema','h','h^2')
