If you don't want to keep around the result of each step, then your algorithm is simply:
A = [1 2 3 4; 5 6 7 8 ; 9 10 11 12; 13 14 15 16];
u = ones (4,1);
for s = 1:10
z = A * u;
m = max(z);
u = z / m;
end
If you do want to keep the intermediary results, then you can indeed use cell arrays, like this:
A = [1 2 3 4; 5 6 7 8 ; 9 10 11 12; 13 14 15 16];
uinitial = ones (4,1);
numsteps = 10;
u = cell(1, numsteps+1); z = cell(1, numsteps); m = cell(1, numsteps);
u{1} = uinitial;
for s = 1:numsteps
z{s} = A * u{s};
m{s} = max(z{s});
u{s+1} = z{s} / m{s};
end
Or you could use arrays which will be easier to inspect
A = [1 2 3 4; 5 6 7 8 ; 9 10 11 12; 13 14 15 16];
uinitial = ones (4,1);
numsteps = 10;
u = zeros(4, numsteps+1); z = zeros(1, numsteps); m = zeros(1, numsteps);
u(:, 1) = uinitial;
for s = 1:numsteps
z(s) = A * u(:, s);
m(s) = max(z(s));
u(:, s+1) = z(s) / m(s);
end
