You have to specify the order. One thing you can do is to specify a high order and return the reflection coefficients. The negative of the reflection coefficients is the partial autocorrelation function. The partial autocorrelation function will decay to zero at what is essentially the optimal AR order.
For example, I create an AR(2) process, but fit an AR(15) model. Then I'll look at the partial autocorrelation function.
A = [1 1.5 0.75];
x = filter(1,A,randn(1000,1));
[arcoefs,E,K] = aryule(x,15);
pacf = -K;
lag = 1:15;
stem(lag,pacf)
You see that the PACF essentially decays to zero after lag 2. That shows that an AR(2) model would be the best.
If you want you can easily construct approximate 95% confidence intervals to help in determining when the PACF values are not significantly different from 0.
stem(lag,pacf,'markerfacecolor',[0 0 1]);
xlabel('Lag'); ylabel('Partial Autocorrelation');
set(gca,'xtick',1:1:15)
lconf = -1.96/sqrt(1000)*ones(length(lag),1);
uconf = 1.96/sqrt(1000)*ones(length(lag),1);
hold on;
line(lag,lconf,'color',[1 0 0]);
line(lag,uconf,'color',[1 0 0]);
