Correct multiple polynomial regression formula with OLS in Python

Question

I need to perform multiple polynomial regression and obtain statistics, p value, AIC etc. As far as I understood I can do that with OLS, however I found only a way to produce a formula using one independent variable, like this:

model = 'act_hours ~ h_hours + I(h_hours**2)'
hours_model = smf.ols(formula = model, data = df)

I tried to define a formula using two independent variable, however I could not understand if that is the correct way and if the results are reasonable. The line that I doubt is model = 'Height ~ Diamet + I(Diamet**2) + area + I(area**2). The full code is this one:

import pandas as pd
import statsmodels.formula.api as smf

train = pd.read_csv(r'W:\...file.csv')

model = 'Height ~ Diamet + I(Diamet**2) + area + I(area**2)'
hours_model = smf.ols(formula = model, data = train).fit()
print(hours_model.summary())

The summary of the regression is here:

                            OLS Regression Results                            
==============================================================================
Dep. Variable:                 Height   R-squared:                       0.611
Model:                            OLS   Adj. R-squared:                  0.609
Method:                 Least Squares   F-statistic:                     376.0
Date:                Fri, 04 Feb 2022   Prob (F-statistic):          1.33e-194
Time:                        08:50:17   Log-Likelihood:                -5114.6
No. Observations:                 963   AIC:                         1.024e+04
Df Residuals:                     958   BIC:                         1.026e+04
Df Model:                           4                                         
Covariance Type:            nonrobust                                         
==================================================================================
                     coef    std err          t      P>|t|      [0.025      0.975]
----------------------------------------------------------------------------------
Intercept         13.9287     60.951      0.229      0.819    -105.684     133.542
Diamet             0.6027      0.340      1.770      0.077      -0.066       1.271
I(Diamet ** 2)     0.0004      0.002      0.262      0.794      -0.003       0.004
area                3.3553      5.307      0.632      0.527      -7.060      13.771
I(area** 2)        0.2519      0.108      2.324      0.020       0.039       0.465
==============================================================================
Omnibus:                       60.996   Durbin-Watson:                   1.889
Prob(Omnibus):                  0.000   Jarque-Bera (JB):               86.039
Skew:                           0.528   Prob(JB):                     2.07e-19
Kurtosis:                       4.015   Cond. No.                     4.45e+05
==============================================================================
Notes:
[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.
[2] The condition number is large, 4.45e+05. This might indicate that there are
strong multicollinearity or other numerical problems.