Inverse Fermi Dirac Integral using Octave/Matlab

I need to obtain the inverse of the fermi dirac integral of order -1/2 in my code. This is what I do. I create a file inv_fermi_dirac_mhalf.m

function y=inv_fermi_dirac_mhalf(fval,init_guess)
global fvalfermi;
%get length of vector
Nel=length(fval);

%initialize
y=zeros(Nel,1);

%loop over all data
for nloop=1:Nel
fvalfermi=fval(nloop);
y(nloop)=fzero('myfunfermi',init_guess(nloop));
end
end

The loop enable the function to take vector data. It uses optimization function fzero to obtain the inverse of the fermi dirac integral. You will need to install Optim and Miscellaneous package into your octave before using fzero function. This function calls another function which will find its zero. So in this case myfunfermi.m is

function y=myfunfermi(x)
global fvalfermi;

y=fermi_dirac_mhalf(x)-fvalfermi;

Note that I need to use a global function to pass into myfunfermi.m. This is because the function that fzero can call must be in the form

And the function fermi_dirac_mhalf.m is from GSL package of Octave.