You can use fmincon,
fmincon(@(x) -x(3), [Y0,L0,z0], [],[],[],[],lb,ub, @(x) nonlcon(x,P, W, iota , nu));
function [c,ceq] = nonlcon(x,P, W, iota , nu)
[Y,L,z]=deal(x(1),x(2),x(3));
ceq = z - P*Y - W*L - z^(-iota)*L^(nu); c=[];
end
but with such a nonlinear constraint, you will need good initial guesses of the unknown parameters Y0,L0,z0 in order to avoid local minima. Because you only have 3 unknowns, it shouldn't be hard to vectorize a discrete grid search for an approximate solution, given appropriate lower and upper bounds lb and ub on the unknowns.
[YY,LL,zz]=ndgrid( linspace(lb(1):ub(1),300) ,...
linspace(lb(2):ub(2),300) , ...
linspace(lb(3):ub(3),300));
con=abs( zz - P.*YY - W.*LL - zz.^(-iota).*LL.^(nu) )<=tolerance;
i=find( zz(con)==max(zz(con)) ,1);
[Y0,L0,z0]=deal(YY(i), LL(i),zz(i)); %approximate solution - intial guess
