What is limiting about a linear model? I've read that a linear model means linear in the parameters, and not necessarily in the predictors. For example, both:
$$Y=\beta_0+\beta_1x_1+\cdots+\beta_kx_k+\epsilon$$
and
$$Y=\beta_0+\beta_1x_1+\beta_2 x_2^2+\beta_3 e^{5x_3}+\cdots+\epsilon$$
are linear models.
Visualy, I would expect such flexibility would let me model any sort of shape between the response and the predictors if I plot my data. I haven't yet learnt more advanced models, but what would be a drawback/incapability of just a linear model like this?
(I realise you wouldn't be able to use linear regression on $Y=\beta_0 + \beta_1 x^{\beta_2}+\epsilon$, for instance, but I'm having trouble visualising/understanding how that would be preventive/inflexible in modelling)
Thanks in advance!
 A: I will cite an educational reference to indicate the possible drawbacks. To quote for the case of a Simple Linear Regression Model:

Objective: model the expected value of a continuous variable, Y, as a linear function of the continuous predictor, X, E(Yi) = β0 + β1xi
Model structure: ${Y_i = β_0 + β_1x_i + \epsilon_i}$
Model assumptions: Y is normally distributed, errors are normally distributed, ${\epsilon_i}$ ∼ N(0, ${σ^2}$), and independent.

In the corresponding case of Generalized Linear Models (GLMs) the assumptions cited include, to quote from the same reference:

The data Y1, Y2, ..., Yn are independently distributed, i.e., cases are independent.
The dependent variable Yi does NOT need to be normally distributed, but it typically assumes a distribution from an exponential family (e.g. binomial, Poisson, multinomial, normal,...)
GLM does NOT assume a linear relationship between the dependent variable and the independent variables, but it does assume linear relationship between the transformed response in terms of the link function and the explanatory variables; e.g., for binary logistic regression ${logit(π) = β_0 + β_X}$.
Independent (explanatory) variables can be even the power terms or some other nonlinear transformations of the original independent variables.
The homogeneity of variance does NOT need to be satisfied. In fact, it is not even possible in many cases given the model structure, and overdispersion (when the observed variance is larger than what the model assumes) maybe present.
Errors need to be independent but NOT normally distributed.
It uses maximum likelihood estimation (MLE) rather than ordinary least squares (OLS) to estimate the parameters, and thus relies on large-sample approximations.

So, the differences from Simple Linear Regression essentially relate to a normality assumption for Y and the error terms, while GLMs do NOT require such an assumption, but do generally operate within the exponential family of distributions.
Also, homogeneity of variance is only in place for Simple Linear Regressions and GLM can specify an appropriate variance-covariance matrix structure.
Lastly, GLMs generally employs a numerically more complex maximum likelihood estimation routine which is not required for ordinary regression.
To answer the particular question: "but what would be a drawback/incapability of just a linear model like this?", the answer is the correct specification of the error structure, and even the diagonal matrix relating to variances, with some explanatory variables involving powers.
A: 
$$Y=\beta_0+\beta_1x_1+\beta_2 x_2^2+\beta_3 e^{5x_3}+\cdots+\epsilon$$
are linear models.
Visualy, I would expect such flexibility would let me model any sort of shape between the response and the predictors if I plot my data. I haven't yet learnt more advanced models, but what would be a drawback/incapability of just a linear model like this?

Yes, you can model any sort of shape.
But the flexibility of the model, as function of the parameters $\beta_i$ is limited. The model parameters only occur in the linear part. So you can't for instance fit this model
$$Y=\beta_0+\beta_1x_1+\beta_2 x_2^{\beta_4} +\beta_3 e^{\beta_5 x_3}+\cdots+\epsilon$$
You can change your model 'shape' $\beta_2 x_2^2+\beta_3 e^{5x_3}$ by changing those coefficients $2$ and $5$ but they are not free model parameters that can be changed in the fitting procedure.

(I realise you wouldn't be able to use linear regression on $Y=\beta_0 + \beta_1 x^{\beta_2}+\epsilon$, for instance, but I'm having trouble visualising/understanding how that would be preventive/inflexible in modelling)

This is a bit of a loaded question. There is not really anything to understand visually. You can make any shape of curve with a linear regression. But multiple shapes will not be available within a single model. For instance you can have the shapes:
$$Y=\beta_0 + \beta_1 x^2+\epsilon$$
or
$$Y=\beta_0 + \beta_1 x^3+\epsilon$$
or using whatever other coefficient.
But only with a more general non-linear model can you capture all those possible shapes at once.
$$Y=\beta_0 + \beta_1 x^{\beta_2}+\epsilon$$
This is for instance useful when the coefficient $\beta_2$ is an unknown parameter that you wish to determine using inference.
A: Just an example: Step functions cannot be represented by linear regressions: A factory at the sea side has a wall to protect it from the waves. Waves smaller then 5 meters stay behind the wall and do no harm. Waves above 5 meters lead to water coming into the cooler, short-circuit it and there is a loss worth of 10 Million Dollars. Model the loss as a function of wave height. 
Simplest problem imaginable for a decision tree regression, not at all a good match for a linear model (even logistic regression claims perfect separation...). 
A: There's little limiting about a linear model per se. In fact, there is Cybenko's universal approximation theorem for neural networks! This has the output of a single layer network as a linear function of some constructed predictors. The problem lies in finding the right set of predictors, generalization out of sample and so forth. In practice these are hard problems.
