% Az egységnégyzeten oldja meg a program a 
% -Laplace u =-20 Poisson-egyenletet
% y=0 esetén u=1-x, x=0 esetén u'_x=y^2-1, 
% különben homogén Dirichlet-peremfeltétellel. 

f=@(x,y)-20+0*(x+y);

% A merevségi mátrix konstrukciója
n=15;
h=1/(n+1);
A=kron(gallery('tridiag',n+1),speye(n))+kron(speye(n+1),gallery('tridiag',n));
A(1:n,1:n)=A(1:n,1:n)/2;

% A terhelési vektor
[xi,yi]=meshgrid(0:h:1-h,1-h:-h:h);
b=h^2*reshape(f(xi,yi),n^2+n,1);
b(1:n)=b(1:n)/2;
b(1:n)=b(1:n)+h*(1-[1-h:-h:h].^2)';
b(2*n:n:n+n^2)=b(2*n:n:n+n^2)+(1-[h:h:1-h])';
b(n)=b(n)+1/2;

% Az egyenletrendszer megoldása 
u=A\b;

% A megoldás ábrázolása
z=zeros(n+2);
z(2:n+1,1:n+1)=reshape(u,n,n+1);
z(n+2,1:n+2)=1-[0:h:1];
tri=[0 0 0];
for i=1:(n+1)*(n+2)
    if mod(i,n+2)~=0
        tri(length(tri(:,1))+1,:)=[i,i+1,i+n+3];
        tri(length(tri(:,1))+1,:)=[i,i+n+2,i+n+3];
    end
end
tri(1,:)=[];
[xi,yi]=meshgrid(0:h:1,1:-h:0);
trisurf(tri,xi,yi,z)
axis equal
xlabel('x','FontSize',14)
ylabel('y','FontSize',14)
title(['n=',num2str(n)],'FontSize',14)