function amath3132hw3
  clear all; close all;
disp('Answer to Q(1)')
disp('==============')
disp('k     (i,j)')
disp(' ')
for i=1:5,
  for j=1:5,
   k=point(i,j,5);
   s=sprintf('%d     (%d,%d)',k,i,j);
   disp(s)
  end
end
disp('      ')
disp('Hit ENTER to continue.....')
pause

disp('Answer to Q(3)')
disp('==============')


nu=5e-1;
n=zeros(4,1);
res=n;
for i=1:length(n)
  n(i)=2^(i+2);
  res(i)=cal_norm(n(i),nu);
end
N=n.^2;
figure(1); clf
loglog(N,res,'o-');
xlabel('N')
ylabel('||b-Ax||_2')


%calculate norm of b-Ax
function [res]=cal_norm(n,nu)
  
  x1=-2;y1=-2;
  xn=2;yn=2;
  
  x=linspace(x1,xn,n);
  y=linspace(y1,yn,n);
  
  dx=max(diff(x));
  dy=max(diff(y));
  
  [x,y]=meshgrid(x,y);
  
  u=exp(-(x.^2+y.^2)./nu);

  [ux,uy]=gradient(u,dx,dy);
  f=divergence(x,y,ux,uy);
  
  for i=1:n
    for j=1:n
      k=point(i,j,n);
      b(k)=f(i,j);
    end
  end
  b=-b.*(dx^2);
  b=b';

  res=norm(b-Au(u,n));

  
%Find Ax, given nxn array u
function [res]=Au(u,n)
% u is a vector
  
  res=zeros(n*n,1);
  
  
  for i=1:n
    if(i==1)
      for j=2:(n-1)
        k=point(i,j,n);
        res(k)=4*u(point(i,j,n))-u(point(i+1,j,n))-u(point(i,j-1,n))-u(point(i,j+1,n));
      end
    elseif(i==n)
      for j=2:(n-1)
        k=point(i,j,n);
        res(k)=4*u(point(i,j,n))-u(point(i-1,j,n))-u(point(i,j-1,n))-u(point(i,j+1,n));
      end
    else %i=2:n-1
      for j=1:n
        k=point(i,j,n);
        if(j==1)
          res(k)=4*u(point(i,j,n))-u(point(i-1,j,n))-u(point(i+1,j,n))-u(point(i,j+1,n));
        elseif (j==n)
          res(k)=4*u(point(i,j,n))-u(point(i-1,j,n))-u(point(i+1,j,n))-u(point(i,j-1,n)); 
        else %j=2:n-1
          res(k)=4*u(point(i,j,n))-u(point(i-1,j,n))-u(point(i+1,j,n))-u(point(i,j-1,n))-u(point(i,j+1,n));
        end    
      end
    end
  end

  
%Q1(a)
function k = point(i,j,n)
  
  k=j+n*(i-1);
  
  
  

