As you have only a limited number of values for f, you should avoid using complicated 3D or 4D arrays (indexing them is not always faster than looping indexing based on simpler 2D arrays), or solutions based on ARRAYFUN, and go for a much simpler solution based on a cell array for storing Y's.
Also, your matrices are quite simple so you could define them in a simpler way that you are doing above, i.e. you seem to have (for any given value of f):
A1 = f + [0 1; 2 3]
A2 = A1 + 1
B = f + [6; 7]
Using a cell array for storing results Y, you could implement the following:
f = [20 500 10000] ;
n = numel(f) ;
Y = cell(n,1) ;
for k = 1 : n
A1 = f(k) + [0 1; 2 3] ;
Y{k} = A1 * (A1 + 1) * (f(k) + [6; 7]) ;
end
Solutions could be then accessed as follows:
Y{1}
ans =
48966
53738
Y[2}
...
etc
Note that this solutions is flexible, as the size of the cell array and the upper bound of the loop adapt to the number of elements in the f vector. Also, I assumed that you just need Y's in the end and not successive values of A1, A2, and B. You could use a few extra cell arrays if you needed to store them.
