There is an obvious singularity at x == 0. This suggests that a polynomial model will be useless, or at best poor. Polynomials do not have singularities, so you would need a high order polynomial. (When you try to fit a polynomial to a function with a singularity, it will do strange things. You will need to use a high order, and that in turn will be a problem.)
Likewise, my own SLM will have problems. Again, splines don't like singularities.
However, IF you swap the x and y axes, you will be able to gain a very nice fit, using a variety of tools. Essentially you will build a model of the form
This model will have no singularities. A spline model (like SLM) will be trivial to fit. Polynomials will also work well enough in this inverse form since there will be no singularity.
To evaluate the model now is slightly harder, but still easy enough. For a polynomial mode, f(y), any given x to find y, use roots or fzero. Thus you would find the roots of the polynomial in y
If f(y) was built using SLM, then use SLMEVAL, which can invert a spline model.
If you insist on trying to fit a model of the form
then you will need terms in the model that have a singularity in them. So a logical starting model might have terms like this:
y = a0 + a1/x + a2/x^2 + a2*x + ...
so a polynomial with some terms with negative exponents. My own polyfitn (also found on the file exchange) would work. Regardless, that model may be difficult to fit well using least squares because of the singularity.