theta(delta) = r1*sin(delta)*u; % <-- Added the (delta) indexing
You will also need to change some of your matrix operators to element-wise operators. E.g.,
E = Psph./((r1.^2)*(r2.^2)*4*pi*(u.^2)*sin(delta).*(cos(theta)).^2);
But note that this entire loop:
for delta = [29 30 31 59 60 61 89 90 91]
theta(delta) = r1*sin(delta)*u;
end
can be replaced with these lines:
delta = [29 30 31 59 60 61 89 90 91]
theta = r1*sin(delta)*u;
