Is there a black box that can extract polynomial (quadratic) relationships? I have an equation of the form $f(x)=\alpha x + \beta x^2$ that I want to match to experimental data.
My current approach is to plot the data on a log-log scale and find the slope unity region (by 'eye-ballin') and use this to extract $\alpha$, then fit the rest of the data to extract $\beta$. 
This doesn't feel very robust. Can anybody point me to a procedure, algorithm, or strategy to extract power laws like this? Perhaps something that can be applied to thousands of data points automatically?
 A: The most obvious approach to fitting the function you specified (quadratic through the origin) would be linear regression. (As I mentioned in comments I would not call that a power law.)
You can specify a model which has $E[Y|x]=\alpha x + \beta x^2$:
Here are some simulated example data (1000 points, from the model, with Gaussian errors):

Here's fitting a quadratic function through the origin (in R):
> quadfit=lm(y~0+x+I(x^2))
> summary(quadfit)

Call:
lm(formula = y ~ 0 + x + I(x^2))

Residuals:
     Min       1Q   Median       3Q      Max 
-1.69152 -0.32394 -0.00318  0.31763  2.43674 

Coefficients:
       Estimate Std. Error t value Pr(>|t|)    
x      0.343774   0.013270   25.91   <2e-16 ***
I(x^2) 0.177502   0.002235   79.40   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.5024 on 998 degrees of freedom
Multiple R-squared:  0.9945,    Adjusted R-squared:  0.9945 
F-statistic: 8.962e+04 on 2 and 998 DF,  p-value: < 2.2e-16

So the estimates are $\hat{\alpha}=0.3438$ and  $\hat{\beta}=0.1775$. 
Here's the fit:

Now that particular fit (in particular, the calculation of the standard errors) assume a constant variance and independent observations. If you don't have spread that's close to constant (in many situations spread is larger when the mean is larger, especially if observations are restricted to be positive), a different fitting approach would make more sense.
Indeed, if you have dependence, or noise in the $x$ (so the actual $x$ you record is not quite identical to the $x$ that passed through to the response variable) it would also not be the best choice.
The more information you have about how the noise operates, the better (in terms of choosing good ways to estimate the parameters).

Questions from comments:

I'm curious if you would know the name of the underlying algorithm that is used to perform the fit? 

The algorithm that was used to perform the fit here is called QR decomposition, but it's one of many ways to perform linear least squares fits; it is one suitable to a wide range of situations and a good default algorithm. It's used by most statistics packages.

I think it can be done under least-squares, 

Least squares isn't the algorithm, it's the criterion that the algorithm optimizes. As I pointed out in my answer above, it's not automatically the best choice, just an easy one to estimate.

but I'm worried that it might require a good starting point. Is there something implicit? 

No starting point is required for linear least squares. Nonlinear least squares and GLMs require starting points (for two common examples), but linear least squares has an explicit solution. The algorithms for linear least squares simply calculate that solution in a numerically stable manner.
A: Fitting power laws to empirical data may be quite tricky. I strongly recommend to use approach proposed in A. Clauset, C.R. Shalizi, and M.E.J. Newman, Power-law distributions in empirical data, SIAM Review 51(4), 661-703 (2009). The authors not only develop an MLE-based method, but also provide its implementation in Matlab, R, Python, and C++. 
