Why can't we control for squares and quadratics of predicted values in regression? If my goal were prediction (e.g. of propensity scores), why couldn't I control for higher order terms of the model equation? For example, why not estimate the model and then control for squares and cubes of the predicted values in a second stage and get an even better estimate of the correct predicted values (rather than just a test of functional form)? I know the standard errors would be wrong but couldn't I bootstrap the process?
For prediction, yes you can consider the models
$$
y = \beta_0 + \beta_1 x_1 + \cdots + \beta_p x_p + \gamma_2 \hat{y}^2 + \cdots + \gamma_m \hat{y}^m + error,
$$
where $\hat{y}$ represents the first-step OLS fitted values and $m$ is chosen by something like cross validation.

I haven't seen this approach used before. My personal guess is that this approach is not as useful as other commonly used ones (SVM, splines, GAM, etc.). For example, if $p$ is large (in comparison to the number of observations $n$), the first-step OLS may already be overfitting so including $\hat{y}$ is not practical. (Yes, you can use lasso residuals but that's a different story.) If $p$ is small, nonlinearity can perhaps be better handled by splines or even by simply augmenting the equation with quadratic and cubic terms of the features. Some generalized additive models (GAM) are already there too.

My personal experience is that nonlinearity is not so important (for prediction using economic data). It is usually a lot more important to avoid overfitting nicely. To me, your suggestion seems to be useful in some cases but not in many.

That said, I do not want to dissuade you from pursuing this issue, although there is a (high) chance of ending up with the conclusion that it is not very useful given the availability of other methods. BTW, you would already know it, but just in case, Hastie, Tibshirani and Friedman's book (The Elements of Statistical Learning) is helpful.