I suggest using "Generalised Additive Models". These type of models are linear but can treat wild non-linearity. The idea is - e.g. with regression splines - that a number of linear regressions are "stacked", so that they can jointly account for highly non-linear effects.
Here is a Python implementation: https://pygam.readthedocs.io/en/latest/
When you are bound to linear regression (OLS), you can add polynomials to the regression. In this case you simply generate a new "column" in your data frame, containing e.g. $x^2$. You can add this variable to the regression directly because linear regression is additive:
Example:
$$ y = \beta_0 + \beta_1 x + u $$
...can be augmented by a squared term for $x$...
$$ y = \beta_0 + \beta_1 x + \beta_2 x^2 + u $$
... and this also works for $x^3$ (and so on) or you can take $log()$ etc.
With GAM you don't have to decide on how to model non-linearity. That is the great advantage of GAMs. When you stick to OLS, you need to check if the non-linearity (imposed by you) really helps to improve fit and/or prediction.
GAM are very well explained in "Introduction to Statistical Learning" have a look at Chapter 7. There also is Python code for the Labs in the book.