If you cannot solve algebraically (no closed form solution), then fzero is the best bet. Something like this:
% example y range
y = linspace(ymin,ymax);
% Create empty array of n values of same size
nsol = zeros(size(y));
% Loop y values
for i=1:numel(y)
% grab local y
ylocal = y(i);
% Setup fzero function function
myfunc = @(n) Pu(n)/Pa(n) - ylocal;
% Solve for myfunc=0
nsol(i) = fzero( myfunc, nguess );
end
Note that "fzero" is being adapted from solving for 0, to solving for ylocal by subtracting that term. There are probably faster ways to get "myfunc" defined outside the loop.
