function phi = solver( a_E, a_W, a_U, a_D, a_P, b )
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%  matrix solution for phi_P at points P_ij
%%  a_P = coefficient of phi_P at point P_ij
%%  a_E = coefficient of phi_E  at point east of P_ij,   i.e.  (i,j+1)
%%  a_W = coefficient of phi_W  at point west of P_ij,   i.e.  (i,j-1)
%%  a_U = coefficient of phi_U  at point up from P_ij,   i.e.  (i-1,j)
%%  a_D = coefficient of phi_D  at point down from P_ij, i.e.  (i+1,j)
%%
%%  system of equations to be solved for phi
%%
%%       big_A * phi = rhs
%%
%%  big_A = sparse matrix with diagonals -a_P, a_E etc, 
%%  phi   = unknown solution at grid points P
%%  rhs   = column vector of constant terms  derived from b
%%             (includes constant source, BC, previous step, etc)
%%
%% On entry:
%%  a_E, a_W, a_U, a_D, a_P and rhs_stuff are matrices with  
%%  the same dimensions as solution matrix.
%%  They have already been modified to include Boundary Conditions
%%  Each entry X_ij in one of these matrices "X" contains terms 
%%  associated with control-volume center point P_ij.
%%
%%  a_P phi_P = a_E phi_E + a_W phi_W + a_U phi_U + a_D phi_D + b
%%
%% or:
%%
%%  - a_P phi_P + a_E phi_E + a_W phi_W + a_U phi_U + a_D phi_D = - b
%%
%%
%%  Each matrix a_E, a_W, a_U, a_D, a_P is reshaped into a
%%  column vector, which becomes a diagonal of the
%%  sparse matrix big_A
%%    a_P         -> main diagonal
%%    a_U and a_D -> adjacent diagonals
%%    a_E and a_W -> diagonals nz+1 above and below main diagonal 
%%
%%  matrix b contains rhs stuff 
%%     - source contributions S_C independent of phi
%%     - B.C. contributions
%%     - values of solution at previous step (if transient problem)
%%  b is reshaped into a column vector rhs
%%         rhs  = right-hand-side vector
%%
%%-----------------------------------------------------------------
%  find nz, nx (number of rows and columns)
         [nz, nx] = size( b );
         n_volumes = nx*nz;
%
% reshape a_X matrices into diagonals of big_A matrix
%   column 3 is center diagonal of big_A
           Diags = zeros( n_volumes, 5 );
       %%    cols = [ -(nz+1) -1  0  1  (nz+1) ];
           cols = [ -nz -1  0  1  nz ];
%
           Diags(:, 5) = reshape(a_W, n_volumes, 1 );
           Diags(:, 4) = reshape(a_U, n_volumes, 1 );
           Diags(:, 3) = reshape(-a_P, n_volumes, 1 ); % note minus sign
           Diags(:, 2) = reshape(a_D, n_volumes, 1 );
           Diags(:, 1) = reshape(a_E, n_volumes, 1 );
%
%   note that spdiags stuffs the transpose of big_A, because
%   all diagonal elements in column 1 refer to grid point (1,1)
%   and I want column 1 to transpose to the row equation for phi(1,1)
           big_A = spdiags( Diags, cols, n_volumes, n_volumes)';
%
%                
% reshape rhs matrix into column vector
         rhs = reshape(-b, n_volumes, 1);   % note minus sign on b
%
%
%        Solve for phi
          phi = big_A \ rhs;
%
%        Reshape column vector into matrix
          phi = reshape( phi, nz, nx );
%    
%
%