How to get a fine grid from coarse grid containing boundaries ?

I am modelling arbitrary layered subsurface properties. I created a matrix U[m,n,l] in which each value determines in which layer is contained. For example M(1,1,1)=1, means that the first cell is in group 1. I have up to 5 groups. In addition, I created 3 arrays X[m], Y[n], Z[l] containing the position for all the points in space. The distance between points is 50 in x direction, 60 in y direction and 70 in

I need now to extract part of the Matrix and do refinement. I created 3 new arrays I[p],J[q] and K[r]. The range of values from these new arrays are within the X,Y,Z array, respectively.

The distance between the position in the I array grid is 10, in J is 12 and in K is 14. Now I need to calculate a new matrix S[p,q,r] in which is hold to which group these points belong.

I created the next function, but it is not correct properly because it holds the old boundary and not a more refined one. In this function for each point I take the closest value from the U Matrix.

I am trying to explain more or less why these function doesn't solve the task correctly in a picture.

if true
function S = interpunit(X,Y,Z,U,I,J,K)
[xnewmesh,ynewmesh,znewmesh]=ndgrid(xnew,ynew,znew);
  vxyz(size(xnewmesh,1)*size(ynewmesh,1)*size(znewmesh,1),1)=99999;
  for k=1:size(znewmesh,1)
      for j=1:size(ynewmesh,1)
          for i=1:size(xnewmesh,1)

                for l=1:size(X,1)

                       if (X(l)-xnewmesh(i))>0
                            mx2=l;
                            mx1=mx2-1;
                            break
                        end
                end

                for m=1:size(Y,1)
                    if (Y(m)-ynewmesh(j))>0
                        my2=m;
                        my1=my2-1;
                        break
                    end
                end

                for n=1:size(Z,1)
                    if (Z(n)-znewmesh(k))>0
                        mz2=n;
                        mz1=mz2-1;

                        break
                    end
                end

    x1=X(mx1);
    xp=xnewmesh(i);
    x2=X(mx2);
    y1=Y(my1);
    yp=ynewmesh(j);
    y2=Y(my2);
    z1=Z(mz1);
    zp=znewmesh(k);
    z2=Z(mz2);

    v000=U(mx1,my1,mz1);
    v100=U(mx2,my1,mz1);
    v010=U(mx1,my2,mz1);
    v110=U(mx2,my2,mz1);
    v001=U(mx1,my1,mz2);
    v101=U(mx2,my1,mz2);
    v011=U(mx1,my2,mz2);
    v111=U(mx2,my2,mz2);

prop=[(x2-xp)*(y2-yp)*(z2-zp) (xp-x1)*(y2-yp)*(z2-zp) (x2-xp)*(yp-y1)*(z2-zp) (xp-x1)*(yp-y1)*(z2-zp) (x2-xp)*(y2-yp)*(zp-z1) (xp-x1)*(y2-yp)*(zp-z1) (x2-xp)*(yp-y1)*(zp-z1) (xp-x1)*(yp-y1)*(zp-z1)];
[a,b]=min(prop);
v=fliplr([v000 v100 v010 v110 v001 v101 v011 v111]);
S(i+(j-1)*size(xnewmesh,1)+(k-1)*size(xnewmesh,1)*size(ynewmesh,1))=v(b);

            end
        end
    end
    S=reshape(S,size(U,1),size(U,2),size(U,3));
end
  end