function u = fish2d(f)% Poisson solver in 2D based on matlab fft% square geometry, homogeneous Dirichlet boundary conditions% requires uniform mesh (2^k-1)^2 interior nodes %% C. T. Kelley, January 6, 1994%% This code comes with no guarantee or warranty of any kind.%% Solves - u_xx - u_yy = f%% Requires sintv.m and isintv.m%% function u = fishsd(f)%%global fish_lal fish_ual ua;n=length(f);m=sqrt(n);hh=.5/(m+1);h2=(m+1)*(m+1);e=ones(m,1);%% Don't compute eigenvalues, factor tridiagonal matrices% or allocate space for the 2-D representation of the solution% more than once. Note the use of global variables to store% fish_ual, fish_lal, and ua.%if exist('fish_ual')==1    [f1, f2]=size(fish_ual);    m1=sqrt(f1);else    m1=m;endif (exist('fish_ual') == 0) | (m1 ~= m)%% allocate room for ua, the 2D representation of the solution%    ua=zeros(m,m);%% set up diagonal of system using eigenvalues of 1D operator%    x=1:m; x=pi*hh*x'; lam=sin(x); lam=h2*(2*ones(m,1) + 4*lam.*lam);    en=h2*ones(n,1); fn=en; gn=en; dx=zeros(n,1);    for i=1:m        fn(i*m)=0;        gn((i-1)*m+1)=0;        dx((i-1)*m+1:i*m)=lam(i)*ones(m,1);    end    al=spdiags([-fn, dx, -gn], [-1,0,1], n,n);%% factor the tridiagonal matrix for solution in the x direction%    [fish_lal,fish_ual]=lu(al);end%% Map f into a rectangular array.% Take sine transform in y direction.%ua(:)=f;ua=sintv(ua).';%% Solve in the x direction%ua(:)=fish_ual\(fish_lal\ua(:));%% Take inverse transform.%ua=isintv((ua.'));%% Map uz into a linear array.%u=ua(:);