3D Plot help

I have a 3D surface plot that looks like this

flux = linspace(0,10,50);
depth = linspace(0,20,50);
[X,Y] = meshgrid(flux,depth);
kd = (-0.001+0.0153).*Y;
a = zeros(1,length(X)) + 2.2921;
density = [2 5 8 12];
for k = 1:4
    kp = -0.1617.*log(density(k));
    beta = -0.0642.*log(density(k))-0.0392;
    SOTE = (a + kd + kp).*X.^beta.*Y;
    figure(4)
    subplot(2,2,k)
    surf(X, Y, SOTE)
    title(['Density is ',num2str(density(k)),'ft']);
    xlabel('Flux (scfm) X');
    ylabel('Depth (ft) Y');
    zlabel('SOTE (%) Z');
    colorbar
    grid on
end

When executed the plot returns this (Don't mind my mis-labeling...)

The problem I'm having is that the x and y axes should have a line of zeros. If the depth is zero, the solution is zero, and the same goes for the flux. I thought the issue could be when meshgrid is applied the X has a column of zeros while the Y has a row of zeros, but I should still have my surface plot hugging the Depth axis. If you just run for one plot and inspect further you can see that the flux axis of the plot doesn't reach zero. Intuition would say that my surface plot would have the front three corners on the floor while the back corner is in the air, and the overall shape would be dome like.

[Jesse Finnell]

The formula has the issue of approaching infinity as flux approaches zero. At zero though, the solution is zero. This discontinuity is the problem that the plot is facing I'm sure.

[darova]

I agree. When you are trying to raise some value to the negative power:

(if b is positive)