close all

clear all

clc

nx = 50;

N = nx*nx;

x = linspace(0,1,nx);

y = x;

[X,Y] = ndgrid(x,y);

dx = x(2)-x(1);

b = -1.36364*X.^4.5 + X.^3.5 - 0.818182*X.^9 + 0.818182*X.^8 - 0.0808081;

b = b * 247.5 .* (Y - 1).*Y.*cos(X);

b = b -45 *(-X.^5.5 + X.^4.5 - 0.444444*X + 0.222222).*(Y - 1).* Y.* sin(X);

b = b -20 *(X - 1).*X.*cos(X) ;

b = exp( X.^4.5 ) .* b;

b = b.*dx.*dx;

isDirichlet = (X==0) | (X==1) | (Y==0) | (Y==1);

b(isDirichlet) = 0.0;

Ad = zeros(N,5);

Ad(:,1) = -cos(X(:));

Ad(:,2) = -cos(X(:)-dx*0.5);

Ad(:,4) = -cos(X(:)+dx*0.5);

Ad(:,5) = Ad(:,1);

Ad(:,3) = -sum( Ad(:,[1,2,4,5]), 2 );

idx = find(isDirichlet(:));

Ad(idx,:) = 0.0;

Ad(idx,3) = 1.0;

A = spdiags(Ad,[-nx,-1,0,1,nx],N,N)';

uh=A\b(:);

uh = reshape(uh,[nx,nx]);

u_exact = 10*X.*Y.*(1-X).*(1-Y).*exp(X.^(4.5));

u_error = abs( (uh - u_exact)./u_exact );

subplot(1,3,1); view(0,90); hold on; colorbar;

surf(x,y,uh,'edgecolor','none');

xlabel('x'); ylabel('y')

title('Computed')

shading interp

axis square

subplot(1,3,2); view(0,90); hold on; colorbar;

surf(x,y,u_exact,'edgecolor','none');

xlabel('x');

title('Exact')

shading interp

axis square

subplot(1,3,3); view(0,90); hold on; colorbar;

surf(x,y,u_error,'edgecolor','none');

xlabel('x');

title('Error')

shading interp

axis square