% Use of FFT for FD approximation to u for
%
%      u_xx + u_yy  = 4,  0 < x,y < 1
%      u(0,y) = y^2, u(1,y) = 1 + y^2, u(x,0) = x^2, u(x,1) = 1 + x^2
%
%  with h = 1/8. This eqn has exact solution uex = x^2 + y^2
%  For simplicity, this script uses FFTs in both directions, rather than 
%  an FFT in one direction and a tridiagonal solver in the other.
%  Note that u_{i,j} corresponds to ROW i and COLUMN j of a matrix of u's;
%  to see the gridpoint values in their correct geometrical orientation
%  so increasing x runs to the right and increasing y upward, 
%  you must display u(:,m:-1:1)'.  Note that in this case, the FD method
%  has zero truncation error since u is parabolic, so the FD soln is exact.
%  Script makes a surface plot of u at the interior grid points.

  L = 1;        % Domain size
  Ne = 16;      % Extended-domain gridpoints...should be a power of 2
  h = 2*L/Ne;   % Grid spacing
  m = Ne/2 - 1; % Number of interior gridpoints in original domain

  x = h*(1:m);
  y = h*(1:m);
  
  f = 4*ones(m,m);
  
  % Correct for Dirichlet BCs
  u0y = y.^2;
  u1y = 1 + y.^2;
  ux0 = x.^2;
  ux1 = 1 + x.^2;
  f(:,1) = f(:,1) - ux0'/h^2;
  f(:,m) = f(:,m) - ux1'/h^2;
  f(1,:) = f(1,:) - u0y/h^2;
  f(m,:) = f(m,:) - u1y/h^2;
  
  g = [zeros(m,1) f zeros(m,1) -f(:,m:-1:1)];    %  Odd extension of f in x
  g = [zeros(1,Ne); g; zeros(1,Ne); -g(m:-1:1,:)]; %  Odd extension of g in y
  ghat = fft2(g);
  [I J]  = meshgrid(0:(Ne-1),0:(Ne-1));
  mu = (4/h^2)*(sin(I*pi/Ne).^2 + sin(J*pi/Ne).^2);
  mu(1,1) = 1;   % Avoid 0/0 division; vhat(1,1) is known a priori to be 0
  v = real(ifft2(-ghat./mu));
  u = v(1:(m+2),1:(m+2)); % Extract u(i,j) from its odd extension
  
%  Plot out solution in interior and print out max-norm error

  [X, Y] = meshgrid(x,y);
  uex = X.^2 + Y.^2;
  err = abs(u(2:(m+1),2:(m+1)) - uex); 
  err = norm(err(:),inf)
  surf(X,Y,u(2:(m+1),2:(m+1)))
  xlabel('x')
  ylabel('y')
  zlabel('u')
  title('FD/FFT for 2D Poisson Eqn')
