Dealing with non-linear variable in multiple linear regression model [duplicate]

Question

I have a multiple linear regression model which should explain the variation store price elasticities using consumer characteristics of the market area surrounding a store. Therefore, my dependent variable are 74 different PE. My independent variables are income (log median income), education level, elderly and household size. The problem I am dealing with is that income and PE show a non-linear (parabolic) relationship in my category (beer).

The values of income range from 9.8 to 11.2. Now my professor gave me two options:

1. Tranform the variable (square)

   demoreg <- lm(pe ~ squaredincome + income + educ + hhlarge + age60, data = demo)
   summary(demoreg)
   

   I feel like there is a mistake but how can I interpret the income coefficient(s)?

   Coefficients:
                 Estimate Std. Error t value Pr(>|t|)    
   (Intercept)  64.588934  34.506386   1.872 0.065735     
   sqincome      0.569562   0.315751   1.804 0.075892     
   income      -12.437625   6.592960  -1.887 0.063696     
   educ          0.019858   0.005524   3.595 0.000626    
   hhlarge      -0.029244   0.013254  -2.207 0.030886    
   age60         0.012861   0.005902   2.179 0.032947
   
   Residual standard error: 0.2441 on 65 degrees of freedom
   Multiple R-squared:  0.5755,    Adjusted R-squared:  0.5428 
   F-statistic: 17.62 on 5 and 65 DF,  p-value: 5.376e-11
   

2. Option: Create two variables one for all income levels that show a negative relationship and one for all income levels that show a positive relationship so I get two coefficients that show how much the PE changes if the median income changes 1%.

   PE = β0 + β1* income_low + β2* income_high
   

   As the low is at 10.6, I created a dummy variable:

   data$dummylow  <- ifelse(data$income <= 10.6,1,0)
   data$dummyhigh <- ifelse(data$income >10.6,1,0) 
   data$low       <- ifelse(data$dummylow ==1, demo_clean$income, NA)  #$
   data$high      <- ifelse(data$dummyhigh ==1, demo_clean$income, NA)
   

   If I include now both variables in the regression equation I get the error code:

   demoreg <- lm(pe ~ low + high + educ + hhlarge + age60, data = demo)
   lm.fit(x, y, offset = offset, singular.ok = singular.ok, ...) All cases NA
   

   I want to include both variables in the equation as I need two coefficients, but I don't know how to do it.