How to vectorize extraction of relevant sections from multidimensional matrix?

The purpose of the attached code is to interpolate a function over multiple nodes and assign weights to the interpolated values. As you can see I am using many for loops. I have used reshape and repmat methods for other sections of my script but I am stuck in the lines that I am attaching to this question. I would appreciate if you can provide me some guidance to vectorize this. If the code is messy or the description is not clear please let me know. Thank you.

 nodes = 5;
age_max = [15 1 1];
% Create a gridded interpolant function
F = griddedInterpolant({agp,tcp,pcp,1:max(age_max),1:lu_bdry,1:no_parcels},Vinit,'spline');
% Interpolate at different quadrature nodes (vectors cx_agp, cx_tcp...)
V_hat = F({cx_agp, cx_tcp, cx_pcp, 1:max(age_max), 1:lu_bdry, 1:no_parcels});
% Assign weights to each interpolated values
w_V_hat = zeros*V_hat;
for j = 1:size(cw_agp,1)
   for jj = 1:size(cw_tcp,1)
     for jjj = 1:size(cw_pcp,1)
       w_V_hat(j,jj,jjj,:,:) = V_hat(j,jj,jjj,:,:)*cw_agp(j)*cw_tcp(jj)*cw_pcp(jjj);
     end
   end
end

Answer 1

[Andrei Bobrov]

[a1,a2,a3] = ndgrid(cw_agp,cw_tcp,cw_pcp);
z = a1.*a2.*a3;
s = size(V_hat);
w_V_hat = V_hat.*repmat(z,[1 1 1 s(4:5)]);

or

[a1,a2,a3] = ndgrid(cw_agp,cw_tcp,cw_pcp);
z = a1.*a2.*a3;
w_V_hat = bsxfun(@times,V,z)

[Raymundo Marcos-Martinez]

Thank you very much for taking the time to answer my question. I had figure out the use of ngrid to simplify the code this way

[a1,a2,a3] = ...
ndgrid(cw_agp,cw_tcp,cw_pcp,1:max(age_max),1:lu_bdry,1:no_parcels);
z = a1.*a2.*a3;
w_V_hat = V_hat.*z

which takes practically the same time to run as the one using bsxfun.

Given the size of the arrays that I am using, the option based on repmat consumes a lot of memory and slows down the computation significantly.

Thank you very much for your guidance.