function ax = atx_sq ( x )

%% ATX_sq computes A*x where A is defined by the function call A(i,j,m)
%
%  Modified:
%
%    7 October 2009
%
%  Author:
%
%    William Martin
%
%  Parameters:
%
%    Input, real x(n,1), the column vector to be multiplied.
%
%    Output, real ax(n,1), the product A*x.
%
  n = length ( x );
  m = sqrt(n);
  h = 1.0/(m+1);
  h2 = 1/h^2;
  ax1=zeros(m,m);
  u=zeros(m,m);
  u(:)=x;
  
  for i=1:m
      for j=1:m
          %get the four diagonal bits
          ax1(i,j)=h2/2.0*(4.0*A(i,j,m)*u(i,j));
          
          %get the other bits -- be sure to handle the edges correctly
          if j~=1
              ax1(i,j) = ax1(i,j) - h2/2.0*((A(i,j,m)+A(i,j-1,m))*u(i,j-1) - A(i,j-1,m)*u(i,j));
          else
              ax1(i,j) = ax1(i,j) + h2/2.0*A(i,j-1,m)*u(i,j);
          end
          
          if i~=1
              ax1(i,j) = ax1(i,j) - h2/2.0*((A(i,j,m)+A(i-1,j,m))*u(i-1,j) - A(i-1,j,m)*u(i,j));
          else
              ax1(i,j) = ax1(i,j) + h2/2.0*A(i-1,j,m)*u(i,j);
          end
          
          if j~=m
              ax1(i,j) = ax1(i,j) - h2/2.0*((A(i,j,m)+A(i,j+1,m))*u(i,j+1) - A(i,j+1,m)*u(i,j));
          else
              ax1(i,j) = ax1(i,j) + h2/2.0*A(i,j+1,m)*u(i,j);
          end
          
          if i~=m
              ax1(i,j) = ax1(i,j) - h2/2.0*((A(i,j,m)+A(i+1,j,m))*u(i+1,j) - A(i+1,j,m)*u(i,j));
          else
              ax1(i,j) = ax1(i,j) + h2/2.0*A(i+1,j,m)*u(i,j);
          end
      end
  end

  ax=ax1(:);