function [A,b]=nma_gen2Ddirch(nx,ny,hx,hy,u) %helper function to generate A,b for solver %nx=9; hx=.125; %ny=9; hy=.125; N= (nx-2)*(ny-2); A= lap2d(ny-2,nx-2); b=zeros(N,1); [X,Y]=meshgrid(0:hx:hx*nx,hy*ny:-hy:0); source=@(X,Y) 1*X; f=source(X,Y); for i=2:ny-1 ii=i-1; for j=2:nx-1 jj=j-1; k=(ii-1)*(nx-2)+jj; if i==2 if j==2 b(k)=f(i,j)+u(i,j-1)+u(i-1,j); else if j==nx-1 b(k)=f(i,j)+u(i,j+1)+u(i-1,j); else b(k)=f(i,j)+u(i-1,j); end end else if i==ny-1 if j==2 b(k)=f(i,j)+u(i,j-1)+u(i+1,j); else if j==nx-1 b(k)=f(i,j)+u(i,j+1)+u(i+1,j); else b(k)=f(i,j)+u(i+1,j); end end else if j==2 b(k)=f(i,j)+u(i,j-1); else if j==nx-1 b(k)=f(i,j)+u(i,j+1); else b(k)=f(i,j); end end end end end end function L2 = lap2d(nx,ny) Lx=lap1d(nx); Ly=lap1d(ny); Ix=speye(nx); Iy=speye(ny); L2=kron(Iy,Lx)+kron(Ly,Ix); %------------------------------------- % % % function L=lap1d(n) e=ones(n,1); L=spdiags([e -2*e e],[-1 0 1],n,n);