Nx = 30; Ny = 30; h = 1/(Nx-1);
A = zeros(Ny^2, Nx^2); b = zeros(Nx^2,1); u = zeros(Nx^2,1);
function r= get_b(i,j,h)
x = i*h; y = j*h;
r = 2*(1-6*x.^2).*y^2.*(1-y.^2) + 2*(1-6*y.^2)*x.^2.*(1-x.^2);
end
for j=1:1:Ny
for i=1:1:Nx
l = i+(j-1)*Nx;
if (i!=1 && j!=1 && i!= Nx && j!= Ny)
A(l,l) = 4;
A(l,l-1) = -1; A(l,l+1) = -1;
A(l,l-Nx) = -1; A(l,l+Nx) = -1;
b(l,1) = h^2*get_b(i,j,h);
else # We are at the boundary
A(l,l) = 1;
end
end
end
# Convert A into sparse CSR
A = sparse(A);
b = sparse(b);
u = A\b;
xi = yi = linspace(0,1,Nx); [Xi,Yi] = meshgrid(xi,yi);
U = reshape(u, [Ny, Nx]);
xa = ya = linspace(0,1,50); [Xa, Ya] = meshgrid(xa, ya);
Ua = (Xa.^2-Xa.^4).*(Ya.^4-Ya.^2);
contour(Xi,Yi,U);
contour(Xa,Ya,Ua); colormap("jet");
daspect([1 1 1]);