axes labels not displaying on 3d subplots

I have created 5 3D subplots which include a for loop. However, the labels of the X and Y axes are not displaying for some reason. I would appreciate any help on this matter. Below is the code.

On a separate note, any suggestions to make the figure more aesthetically pleasing would also be much appreciated.

% parameters
b=0.5;
O=27;
x=1:1:5;

% energies 
e1 = 1:1:100;
e2 = 1:1:100; 

% function
[e1,e2]=meshgrid(e1,e2);

hb=@(x)((O.^2)./factorial(O-x)).*...
    exp(-b.*O.*e2);

hu=@(x)(O.^x).*...
    exp(-b.*O.*e1);

p=@(x)hb(x)./(hb(x)+hu(x));

f=figure('visible','on')
clf(f);

for i=x
    subplot(2,3,i);
    mesh(e1,e2,p(i))
    title(['X = ',int2str(i)], 'FontSize',12);

    % log all axes
    set(gca, 'XScale', 'log');
    set(gca, 'YScale', 'log');
    set(gca, 'ZScale', 'log');

    axis([1 100 1 100 10^-300 1])
    axis square
    grid off 
    set(gca,'FontSize',10)

    xlabel('e1')
    ylabel('e2')
    zlabel('p_{H}')
end

Solution

The issue seems to be something internal to MATLAB with how it is setting the position of the x and y labels when a 3D surface plot is used. This doesn't happen with a basic plot3 plot. If you do a get(get(gca,'Xlabel','Position')), you see that the z coordinate of the label is set to infinity, which I would guess is the problem.

I've come up with a less than ideal workaround, but it seems to accomplish the task:

% parameters
b=0.5;
O=27;
x=1:1:5;

% energies 
e1 = 1:1:100;
e2 = 1:1:100; 

% function
[e1,e2]=meshgrid(e1,e2);

hb=@(x)((O.^2)./factorial(O-x)).*...
    exp(-b.*O.*e2);

hu=@(x)(O.^x).*...
    exp(-b.*O.*e1);

p=@(x)hb(x)./(hb(x)+hu(x));

f=figure('visible','on');
clf(f);

for i=x
    subplot(2,3,i);
    mesh(e1,e2,p(i))
    title(['X = ',int2str(i)], 'FontSize',12);

    % log all axes
    set(gca, 'XScale', 'log');
    set(gca, 'YScale', 'log');
    set(gca, 'ZScale', 'log');

    axis([1 100 1 100 10^-300 1])
    axis square
    grid off 
    set(gca,'FontSize',10)
    xlabel('e1')
    ylabel('e2')
    zlabel('p_{H}')
    set(get(gca,'xlabel'),'Units','Normalized','Position',[0.75 0 0])
    set(get(gca,'ylabel'),'Units','Normalized','Position',[0 0.05 0])
end

You'll probably have to manipulate those position vectors to get the labels exactly where you'd like.

I would also submit a bug report and see what MathWorks says.