function fast_membrane(N)

h = 1/N;  M = (N-1)^2;
b = -h^2*ones(M,1);

A = membrs_fast(N); % use sparse data structure for A
u = A\b;            % backslash operator uses optimized algorithm for sparse matrices
                    % (lu does not)

center_displacement = u( (N/2-1)*(N-1) + N/2 )
surf(reshape(u,N-1,N-1))

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function A = membrs(N)
% find matrix A for square membrane divided into N^2 smaller squares
% use sparse data structure

% For large N this is slow in Matlab because of two nested for loops

n = N-1;
M = n^2;

A = sparse(M,M);  %   initializes sparse data structure 
                  %   (does not use memory for zero entries)
for i=1:n
  for j=1:n
    k = i + n*(j-1);    % number of point in row j, column i
    A(k,k) =  1;   
    if i>1
      A(k,k-1) =  -.25;      % point to the left
    end
    if i<n
      A(k,k+1) =  -.25;      % point to the right
    end
    if j>1
      A(k,k-n) =  -.25;      % point below
    end
    if j<n
      A(k,k+n) =  -.25;      % point above
    end
  end
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function A = membrs_fast(N)
% find matrix A for square membrane divided into N^2 smaller squares
% use sparse data structure

% faster version: use spdiags to avoid double for loop in Matlab

n = N-1;
M = n^2;
e = ones(M,1);
A = spdiags([-.25*e,-.25*e,e,-.25*e,-.25*e], [-n,-1,0,1,n],M,M);
for i=1:n-1   % correct wrong values on sub and superdiagonal
  A(i*n,i*n+1) = 0;  
  A(i*n+1,i*n) = 0;
end







