Feature diffmat rectangular and BCs

diffmat@trunk has been augmented by adding rectangular differentiation matrices, Fourier differentiation matrix, and boundary condition functionality. 

New calling sequences:

```
function D = diffmat(N, varargin)
%DIFFMAT   Differentiation matrix.
%   D = DIFFMAT(N) returns the NxN differentiation matrix associated with the
%   Chebyshev spectral collocation method at second-kind Chebyshev points. 
%
%   D = DIFFMAT(N, K) returns the differentiation matrix of order K. See
%   COLLOC2.DIFFMAT for further details.
%
%   D = DIFFMAT(N, K, DOM) scales the differentiation matrix D to the domain
%   DOM. DOM should be a 1x2 vector.
%
%   D = DIFFMAT(N, K, DOM, DISC) or DIFF(N, K, DISC) returns the differentiation
%   matrix associated with the CHEBDISCRETIZATION DISC.
%
%   D = DIFFMAT(N, 'periodic') returns the N x N first-order Fourier 
%   differentiation matrix 1st-order Fourier diff mat.
%
%   D = DIFFMAT(N, K, 'periodic') returns the N x N Fourier differentiation 
%   matrix of order K.
%
%   D = DIFFMAT(N, K, 'periodic', DOM) scales the Kth-order Fourier 
%   differetiation matrix to the domain DOM.
%
%   D = DIFFMAT(N, K, 'rect') or D = DIFFMAT([N -K], K) returns (N-K) x N 
%   rectangular differentiation matrix of order K which maps from an N-point 
%   Chebyshev grid of second kind to an (N-K)-point Chebyshev grid of first kind.
%
%   D = DIFFMAT([N -K], K, DOM) returns (N-K) x N rectangular differentiation 
%   matrix of order K which is scaled to the domain DOM.
%
%   D = DIFFMAT([N -K], K, DOM, DISC) returns the Kth-order(N-K) x N rectangular 
%   differentiation matrix on domain DOM associated with the CHEBDISCRETIZATION 
%   DISC. When DISC is specified as 'colloc1' or colloc1() or @colloc1, D maps
%   between Chebyshev grids of first kind of size N and (N-K). when DISC is set
%   'colloc2' or colloc2() or @colloc2, D is same as default and maps from a 
%   second-kind Chebyshev grid.
%
%   D = DIFFMAT([M N]) returns an M x N first-order rectangular differentiation 
%   matrix which maps from an N-point Chebyshev grid of second kind to an 
%   M-point Chebyshev grid of first kind.
%   
%   D = DIFFMAT([M N], K) returns an M x N rectangular differentiation matrix of 
%   order K which maps from an N-point Chebyshev grid of second kind to an 
%   M-point Chebyshev grid of first kind.
%
%   D = DIFFMAT([M N], K, DOM) returns the same D but scaled to the domain DOM.
%
%   D = DIFFMAT([M N], K, DOM, DISC) returns a rectangular differentiation 
%   matrix associated with the CHEBDISCRETIZATION DISC. D maps from the N-point
%   Chebyshev grid of first kind when DISC is specified as 'colloc1' or 
%   colloc1() or @colloc1. It is same as default when DISC is 'colloc2' or 
%   colloc2() or @colloc2. 
%
%   D = DIFFMAT(N, K, DOM, DISC, LBC, RBC) or D = DIFFMAT([N -K], K, DOM, DISC, 
%   LBC, RBC) or D = DIFFMAT([N-K N], K, DOM, DISC, LBC, RBC) returns a square
%   differentiation matrix which is deflated by including information about the
%   boundary conditions specified by the cell structures LBC and RBC. When the
%   original differentiation matrix (without boundary condition information) is 
%   square, the boundary conditions are included by row replacement, whereas 
%   this is done by row appending when the original differentiation matrix is 
%   rectangular. The string specifier 'd' (or 'D') and 'n' (or 'N') indicate
%   Dirichlet and Neumann boundary conditions respectively. The character or
%   string 's', 'S', 'sum', 'i', 'I' indicates a side condition of definite 
%   integral i.e. sum(U), where U is the solution to the system formed by D.
%
%   Example 1: D = DIFFMAT(N, K, DOM, DISC, {'d'}, {'n'}) replaces the first and
%   the last row of a square differentiation matrix by Dirichlet and Neumann
%   boundary conditions respectively.
%
%   Example 2: D = DIFFMAT([N -K], K, DOM, DISC, {}, {'n' 's'}) appends two rows 
%   to an (N-K) x N rectangular differentiation matrix to square it up. The 
%   first appended row corresponds to Neumann boundary condition at the right 
%   endpoint while the second appended row corresponds to a side condition 
%   sum(U) = I for a scalar I, where U is the solution to the resulting system.
%   
%   Note that for the case of row appending, the number of the boundary 
%   conditions given by LBC and RBC must total K. Also note that RBC can be 
%   omitted if there are only left boundary conditions.
%
% See also DIFF, COLLOC2.DIFFMAT, CUMSUMMAT.
```

This should be `BARYMAT(Y, X, W, S, R)`.