Apologies if this has already been answered, but I have been unable to find a clear answer/I am not exactly sure what I am looking for. I would like to plot nonlinear predicted probabilities from a multivariate model against a continous variable included in the model at two different levels of an indicator variable. I estimate predicted probabilties of students passing an exam. Here's the model I'm fitting:

Code:
 areg pass i.low_income std_prior_test_score std_prior_test_score_sq std_prior_test_score_cubed other_covs, absorb(fixedeffect)
, where pass is a binary dummy variable =1 if a student passed the exam, low-income is a binary dummy =1 if a student is classified as low-income, std_prior_test_score is the student's standardized test score on a prior exam, std_prior_test_score_sq is that score squared and std_prior_test_score_cubed cubed.

I want to plot the predicted probability of passing against std_prior_test_score at the means of the other covariates. I have played around with the lpoly and margins commands, but have run into issues. Here's what I've done with margins:
Code:
margins i.lowinc_max, at(std_prior_test_score_sq(-2.5(.1)1.5)) atmeans
marginsplot
and obtained the following:
Array


Of course, I'm hoping to obtain the nonlinear relationship with the x-axis variable. Am I able to obtain this nonlinear predictions using the margins command? Or perhaps using a combination of predict and lpoly?

Eventually, I would like to plot the predictions from the following model as well:

Code:
 areg pass i.low_income i.low_income#c.std_prior_test_score i.low_income#c.std_prior_test_score_sq i.low_income#c.std_prior_test_score_cubed other_covs, absorb(fixedeffect)
I appreciate any and all input. I'm sure there exists an easy solution that I am just unfamilar with. Thank you.