First, note this
x = [1 6 9]
x.^-1/3
x.^(-1/3)
(x.^-1)/3
Then, you do not need to store z and w. Here is a more efficient solution:
idx = [1 10 100 3684] ;
N = numel(idx) ;
p = zeros(N,2) ; % pre-allocation
for k = 1:N
w = idx(k) ;
x = R(1:12,w) ;
xz = x.^-1/3;
y = S(1:12,w);
p(k,:) = polyfit(xz,y,1);
f = polyval(p,xz) ;
plot(x, y, 'o', x, f, '-')
end
If you want to store the fitted values as well, pre-allocate f
f = zeros(12,N) ;
...
f(:,k) = polyval(p,xz) ;
