% poissondirichletdf2d.m
% Metodo de diferencias finitas en 2d
% para el problema
%
% - k (u_xx + u_yy) = f(x,y),    (x,y) en Omega = (a,b) x (a,b)
%            u(x,y) = alfa(x,y), (x,y) en la frontera de Omega
%

%% Parametros del problema
a=0; b=1;
k = 0.001;
f = @(x,y) ones(size(x));
alfa = @(x,y) zeros(size(x));

%% Parametros de la resolucion
N = 30;

%% Resolucion
h=1/N;
xx=linspace(a,b,N+1);
yy=xx;
[x,y] = meshgrid(xx,yy);
x = x'; x=x(:);
y = y'; y=y(:);
F=f(x,y);
ALFA=alfa(x,y);

diagonales = ones((N+1)^2,1)*(k*[ -1 -1 4 -1 -1 ]/h^2);

M = spdiags(diagonales, [N+1 1 0 -1 -(N+1)], (N+1)^2, (N+1)^2);
M = M';

rhs = F;

% Ahora corregimos las ecuaciones de los nodos del borde
% piso techo pared izq pared der
ind=[1:N+1 (N+1)^2-N:(N+1)^2  N+2:N+1:N^2  2*(N+1):N+1:(N+1)*N];

M(ind,:) = 0;
rhs(ind) = ALFA(ind);

M = M + sparse(ind,ind,1);
% for k = ind
%     M(k,k)=1;
%     rhs(k)=alfa(k);
% end

U = M\rhs;

U = vec2mat(U,N+1);

mesh(xx,yy,U);
xlabel('x')
ylabel('y')
title(sprintf('N = %d', N))