A = 1000;
n = 12;
r = .1:.01:.2;
k = (15 : 5 :25).';
t1 = bsxfun(@(r,k)(1 + r/n) .^(n*k),r,k);
P1 = A*bsxfun(@times,r,t1 ./ (n*(t1-1)));
P = reshape([100*r;P1],1,[])
variant with one cycle
A = 1000;
n = 12;
r = .1:.01:.2;
k = 15 : 5 :25;
nr = numel(r);
nk = numel(k)+1;
out = zeros(1,nr + (nk-1)*nr);
for ii = 1:nr
t = (1 + r(ii)/n) .^(n*k);
out(nk*ii+(-3:0)) = [100*r(ii), A*r(ii)*t./(n*(t-1))];
end
