%% Poissons ekvation på [0,1]^2
% Detta Matlab-program löser u_xx + u_yy = f(x,y) på [0,1]^2
f = @(x,y) -pi^2*(x.^2+y.^2).*sin(pi*x.*y);

n = 41;     % Antalet punkter i vardera x och y
h = 1/(n-1);

x = 0:h:1;  % Punkter i x och y
y = 0:h:1;

[X,Y] = ndgrid(x,y);    % Skapa repetitioner
xy = [X(:) Y(:)];       % xy(k,:) = [x(i) y(j)], där k = i + (j-1)*n

%% Skapa vektorer som indexerar randen
bottom = 1:n;
top = (n*(n-1)+1):n^2;
left = 1:n:(n*(n-1)+1);
right = n:n:n^2;
% Hela randen ges av alla unika element i ovanstående vektorer
bdry = unique([bottom top left right]);
% Resterande element hamnar i inner (skillnaden mellan alla och bdry)
inner = setdiff(1:n^2,bdry);

%% Skapa finit differensmatris
e = ones(n^2,1);
A = spdiags(1/h^2*[e e -4*e e e],[-n -1 0 1 n],n^2,n^2);

%% Ställ upp ekvationssystem
u = zeros(n^2,1);
% Randvillkor
u(bottom) = 0;
u(top) = sin(pi*x);
u(left) = 0;
u(right) = sin(pi*y);
% Högerledet i ekvationssystemet (bara för inre punkter!)
rhs = f(xy(inner,1),xy(inner,2)) - A(inner,bdry)*u(bdry);

% Löser för de inre punkterna
u(inner) = A(inner,inner)\rhs;

%% Plotta lösningen
mesh(X,Y,reshape(u,size(X)));
xlabel('x'), ylabel('y'), zlabel('u')

%% Skriv ut felet
u_exact = sin(pi*xy(:,1).*xy(:,2));
err = max(abs(u-u_exact));
disp(['h = ' num2str(h) ', err = ' num2str(err) '.'])