
clear all
grid_no  = input('mesh number = ');
% x        = linspace(0, 1, grid_no);
% y        = linspace(0, 1, grid_no);
leng_h   = 1 / ( grid_no + 1 );
x        = leng_h:leng_h:grid_no*leng_h;
y        = x;
[xx, yy] = meshgrid(x, y);
zz       = sin(pi*xx).*sin(pi*yy);

tic;
mtx_A    = mtx_generate_2DLaplacian(grid_no);
rhs_f    = Construct_rhs(grid_no);
sol      = mtx_A \ rhs_f;
fprintf('CPU time for solving linear system = %14.4e, matrix size = %7.0f \n', toc, size(mtx_A,1));

mtx_sol  = zeros(grid_no);
for ii = 1:grid_no
    mtx_sol(ii,:) = sol((ii-1)*grid_no+1:ii*grid_no);
end

figure(1)
mesh(xx, yy, zz);
xlabel('x');
ylabel('y');
zlabel('u');
title('u(x,y) = sin( x \pi ) sin( y \pi )');

figure(2)
surf(xx, yy, zz);
shading interp;
xlabel('x');
ylabel('y');
zlabel('u');
title('u(x,y) = sin( x \pi ) sin( y \pi )');

figure(3)
surf(xx, yy, mtx_sol);
shading interp;
xlabel('x');
ylabel('y');
zlabel('u');
title('Approximate solution');

figure(4)
surf(xx, yy, abs(mtx_sol-zz));
% shading interp;
xlabel('x');
ylabel('y');
zlabel('u');
title('Absolute Error');
