p36.m

% p36.m - Laplace eq. on [-1,1]x[-1,1] with nonzero BC's
% Set up grid and 2D Laplacian, boundary points included:
N = 24; [D,x] = cheb(N); y = x;
[xx,yy] = meshgrid(x,y); xx = xx(:); yy = yy(:);
D2 = D^2; I = eye(N+1); L = kron(I,D2) + kron(D2,I);
% Impose boundary conditions by replacing appropriate rows of L:
b = find(abs(xx)==1 | abs(yy)==1);
% boundary pts
L(b,:) = zeros(4*N,(N+1)^2); L(b,b) = eye(4*N);
rhs = zeros((N+1)^2,1);
rhs(b) = (yy(b)==1).*(xx(b)<0).*sin(pi*xx(b)).^4 + ...
.2*(xx(b)==1).*sin(3*pi*yy(b));
% Solve Laplace equation, reshape to 2D, and plot:
u = L\rhs; uu = reshape(u,N+1,N+1);
[xx,yy] = meshgrid(x,y);
[xxx,yyy] = meshgrid(-1:.04:1,-1:.04:1);
uuu = interp2(xx,yy,uu,xxx,yyy,'cubic');
clf, subplot('position',[.1 .4 .8 .5])
mesh(xxx,yyy,uuu), colormap([0 0 0])
axis([-1 1 -1 1 -.2 1]), view(-20,45)
text(0,.8,.4,sprintf('u(0,0) = %12.10f',uu(N/2+1,N/2+1)))