Try to use a cell array to store the result as follows:
zn = 3; % no. of planes
zmax = 200; % z max
zmin = 100; % z min
zA = linspace(zmax,zmin,zn); % linearly distributed z planes between z max & z min
phiA = (-180:180); % phi value from phi min to phi max
thetaA = (-180:180); % theta value from theta min to theta max
for zi = 1:length(zA)
z = zA(zi);
Pall = []; % Initialize
for phii = 1:length(phiA)
phi = (phiA(phii));
for thetai = 1:length(thetaA)
theta = (thetaA(thetai));
px = phi + theta;
py = phi - theta;
P = [px, py]; % Store value for px and py only
if px>-30 && px<30 && py>-30 && py<30
Pall = [Pall;P]; % Append the point which satisfy the condition
end
end
end
S{zi} = Pall; % Save into a cell array
fprintf('Processing %d of %d ...',zi,zn); % progress
fprintf('done.\n'); % display when done
end
S
