clear A;               // Remoção de valores antigos das variáveis. 
clear b;

stacksize(200000000); // Aumento da memória alocada.

// Dimensões do domínio da equação e pontos da malha para aproximação.
Lx=10; n=101;
Ly=3;  m=101;

hx=Lx/(n-1);
hy=Ly/(m-1);


A=sparse([1 1],[1],[n*m n*m]);

// Termos da matriz associados ás diferenças finitas para pontos internos.
for i=2:n-1
    for j=2:m-1
        A((j-1)*n+i,(j-1)*n+i+1)=hy^2;
        A((j-1)*n+i,(j-1)*n+i-1)=hy^2;
        A((j-1)*n+i,j*n+i)=hx^2;
        A((j-1)*n+i,(j-2)*n+i)=hx^2;
        A((j-1)*n+i,(j-1)*n+i)=-2*(hx^2+hy^2);
        b((j-1)*n+i)=0
    end
end

// Termos da matriz associados ás diferenças finitas para pontos na fronteira.
j=1;
for i=1:n
    A(i,i)=1;
    b(i)=1 + (i-1)*4/(n-1);
end

// Termos da matriz associados ás diferenças finitas para pontos na fronteira.
i=n;
for j=1:m
    A((j-1)*n+i,(j-1)*n+i)=1;
    b((j-1)*n+i)=5 - 2*(j-1)/(m-1);
end

// Termos da matriz associados ás diferenças finitas para pontos na fronteira.
j=m;
for i=n:-1:1
    A((j-1)*n+i,(j-1)*n+i)=1;
    b((j-1)*n+i)=3 - 5*(i-n)/(n-1);
end

// Termos da matriz associados ás diferenças finitas para pontos na fronteira.
i=1;
for j=m:-1:1
    A((j-1)*n+i,(j-1)*n+i)=1;
    b((j-1)*n+i)=8 + 7*(j-m)/(m-1);
end

// Solução do sistema.
z=A\b;


// Reconstrução do solução bidimensional u.
for i=1:n
    for j=1:m
        u(i,j)=z((j-1)*n+i);
    end
end
clear z;

// Montagem do gráfico
fig=scf();
fig.color_map=hotcolormap(32);
grayplot(1:n,1:m,u);
