mnrfit.m

X = [20.5156 12.0474 -12.1288 7.4995 38.9435 -4.5757 20.6246 19.6316 -7.3335 -2.6384 9.1368 24.2808 14.8241 15.9592 5.9249 19.085 -13.4054 -3.3012];
Y = [4       5       5        3      2        3       5       2       3       2      5       4       4       1      6       6       6       6];
[B, dev, stats] = mnrfit(X,Y, 'model', 'ordinal', 'interactions', 'off');

% output:
% B =
%
%    -3.3364
%    -1.6701
%    -0.8469
%    -0.1595
%     0.9495
%     0.0386
%
%
% dev =
%
%    60.5118
% stats
% check  http://www.mathworks.com/help/stats/mnrfit.html#outputarg_stats  for information in stats