The following script fits a quadratic line in the plot (in blue):
quadratic_model <- lm(D ~ GDPCAP + I(GDPCAP^2), data = mydata)
order_id <- order(mydata$GDPCAP)
lines(x = mydata$GDPCAP[order_id],
y = fitted(quadratic_model)[order_id],
col = "blue",
lwd = 2)

Coefficients of the quadratic equation are obtained as follows:
coeftest(quadratic_model, vcov. = vcovHC, type = "HC1")
t test of coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 49.970827 0.310524 160.924 < 2.2e-16 ***
GDPCAP 13.790826 0.340038 40.557 < 2.2e-16 ***
I(GDPCAP^2) -1.180603 0.070727 -16.692 < 2.2e-16 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1