3D interpolation and speed developing predictions

Does anybody have a way to speed this up?

I am doing a 3D interpolation using scatteredInterpolant. The interpolation itself (developing the formula from the sampled data) does not take very long (under 30 seconds). What takes a long time is using the generated formula to develop a matrix of interpolations - 5 minutes.

My data going into the interpolation is around 10,000 x,y locations with a z value. It's as if one walked a field in an evenly spaced grid and took an elevation reading at every point of the grid. I then need to estimate all points that would fill the boxes between these sampled points. Ultimately the final interpolated field data needs to be 55M points.

I am not looping or doing matrix operations on the data at all. The time is all taken up by using the formula to generate the filled in values.

If I could use the sparse matrix of data I would. What I end up doing needs to be nonsparse.

yDiffFunction = scatteredInterpolant(centroids(:,4),centroids(:,3),centroids(:,5));
xDiffFunction = scatteredInterpolant(centroids(:,4),centroids(:,3),centroids(:,6));
xDiff = xDiffFunction(Xq,Yq);
yDiff = yDiffFunction(Xq,Yq);