This is the solution:
X=[wL,wIN1,wIN2,wIN3,wIN1,wIN2,wIN3,wR,wL,wIN1,wIN2,wIN3,wR];
Y=[hIN,hB,hB,hB,hIN,hIN,hIN,hIN,hBC,hBC,hBC,hBC,hBC];
Z=[lL,lB1,lB2,lB3,lIN1,lIN2,lIN3,lR,BC1,BC2,BC3,BC4,BC5];
n=(1:length(Z));
Ztop=(max(Z)+10).*ones(1,length(n));
plot3(app.Map, X(n),Y(n),Ztop,'.k','markersize',12)
view(app.Map, 2)
hold(app.Map, 'on')
[x y] = meshgrid([min(X):1:max(X)],[min(Y):1:max(Y)]);
options = fitoptions('Method','BiharmonicInterpolant');
sf = fit([X',Y'],Z','biharmonicinterp',options);
z = sf(x,y);
surf(x,y,z,'Parent', app.Map,'FaceColor','interp', 'EdgeColor','None')
axis(app.Map, 'equal')
caxis(app.Map, 'auto')
app.Map.YLim = [hB hBC];
app.Map.XLim = [wL wR];
view(app.Map, 2)
colormap(app.Map);
colorbar(app.Map);
hold(app.Map, 'off')