Vectorization of four nested for-loops - bsxfun?

Hi there,

I've been trying to figure out if it is possible to vectorize four nested for-loops of the type shown below. It's a fairly hefty discrete integral computation which takes seriously long time to compute due to the size of the data-set contained in the matrix E. I've seen various threads talking about bsxfun for increasing computational speed, but can't figure out how to do something like this with it, if at all possible?

    parfor l = 1:c1
            for h = 1:c2
                for n = 1:Ny
                    for m = 1:Mx
                        R = sqrt(r0^2+xp(m)^2+yp(n)^2-2*r0*sin(theta(h))*(xp(m)*cos(phi(l))+yp(n)*sin(phi(l))));
                        Er(l,h) = Er(l,h)+sin(theta(h))*sin(phi(l))*E(n,m)*exp(-1i*k*R)/R;
                        Etheta(l,h) = Etheta(l,h)+cos(theta(h))*sin(phi(l))*E(n,m)*exp(-1i*k*R)/R;
                        Ephi(l,h) = Ephi(l,h)+cos(phi(l))*E(n,m)*exp(-1i*k*R)/R;
                    end
                end
                Srad(l,h) = (abs(Etheta(l,h))^2+abs(Ephi(l,h))^2)/(2*Eta);
                Prad = Prad + Srad(l,h)*sin(theta(h))*r0^2*dtheta*dphi;
            end
        end

Any help would be deeply appreciated!

In advance, thank you.

Best Regards,