If you have LSQLIN in the Optimization Toolbox, this can be done with a little bit of effort, as described here http://www.mathworks.com/support/solutions/en/data/1-12BBUC/
% Making some random data
Q1 = 0:0.01:1;
T = cos(2.1*pi*Q1)+0.2*randn(size(Q1));
plot(Q1,T,'k.');
% Polyfit without constraints
order = 4;
C1 = polyfit(Q1,T,order);
hold on;
plot(Q1,polyval(C1,Q1))
V = bsxfun(@power,Q1(:),order:-1:0); % Make Vandermonde Matrix
% Make Constraints on the derivatives
Aleft = [(order:-1:1).*Q1(1).^(order-1:-1:0) 0];
Aright = [(order:-1:1).*Q1(end).^(order-1:-1:0) 0];
Aeq = [Aleft; Aright];
beq = [0;0]; %Value of the derivative is set to zero
% Call LSQLIN with options to prevent warnings
opts = optimset('lsqlin');
opts.LargeScale = 'off';
C2 = lsqlin(V,T,[],[],Aeq,beq,[],[],[],opts);
plot(Q1,polyval(C2,Q1),'r')
hold off;
legend({'Data' 'Polyfit' 'Constrained Polyfit'},'location','best');
