Kelly Hilands

Greg Reeves' version indeed runs much faster, but gave me strange results in tests:
A =
1 2 3
4 5 6
7 8 9
smooth2(A,1) =
1.8 3.6 2.6
4.9 8.9 6.2
4.3 7.6 5.1

I switched rows and columns on the eL and eR matrices and was able to obtain results consistent with the documentation of Kelly Hilands' original code (at least for this test):
smooth2(A,1) =
3 3.5 4
4.5 5 5.5
6 6.5 7

Here is the corrected code (two lines are changed from Greg Reeves' code above)

%Initial error statements and definitions
if nargin<2, error('Not enough input arguments!'), end

N(1) = Nr;
if nargin<3
    N(2) = N(1);
else
    N(2) = Nc;
end

if length(N(1))~=1, error('Nr must be a scalar!'), end
if length(N(2))~=1, error('Nc must be a scalar!'), end

[row,col] = size(matrixIn);
eL = spdiags(ones(row,2*N(1)+1),(-N(1):N(1)),row,row);
eL = eL./(repmat(sum(eL,2),1,row)); % THIS LINE CHANGED FROM REEVES' CODE
eR = spdiags(ones(col,2*N(2)+1),(-N(2):N(2)),col,col);
eR = eR./(repmat(sum(eR,1),col,1)); % THIS LINE CHANGED FROM REEVES' CODE

matrixOut = eL*matrixIn*eR;

I took out the for loops and just used array algebra to do essentially the same thing. Improved speed ~10 fold. I didn't bother checking whether edge behavior is the same. My edge behavior is similar to the default behavior in Matlab's builtin smooth function.

%Initial error statements and definitions
if nargin<2, error('Not enough input arguments!'), end

N(1) = Nr;
if nargin<3
    N(2) = N(1);
else
    N(2) = Nc;
end

if length(N(1))~=1, error('Nr must be a scalar!'), end
if length(N(2))~=1, error('Nc must be a scalar!'), end

[row,col] = size(matrixIn);
eL = spdiags(ones(row,2*N(1)+1),(-N(1):N(1)),row,row);
eL = eL./(repmat(sum(eL,1),row,1));
eR = spdiags(ones(col,2*N(2)+1),(-N(2):N(2)),col,col);
eR = eR./(repmat(sum(eR,2),1,col));

matrixOut = eL*matrixIn*eR;