Weird pdf of a quadratic function of a N(0,1) variable: miscoding or big rounding error? I would like to calculate the pdf of a random variable y defined by :
y=c+b*x+a*x^2
The pdf is a non-central chi-squared distribution. For a>0, it should be equal to zero if y is smaller than d, where d=c-(b^2)/4a.
Strangely enough, when computing it with R, the pdf shoots up at y>d+e, where e is quite large.
Is there an error in my codes (below) or is it a rounding error? In the latter case, how to address it?
set.seed(101)        
x <- seq(-3.5,3.5,length.out=1000)
c<-80   
b<-30   
a<-6    
y<-c+b*x+a*(x^2) # g(x) 
## min(y)

Graph 1: just to get an idea of the function
plot(x[order(x)],y[order(x)],   
type="l",lwd=2, xlim=c(-4,4),   
ylab="y",xlab="x",  
main="a. y=g(x)and density of x")   
par(new=T)  
fx<-exp(-0.5*(x^2))/sqrt(2*pi)  
fx<-dnorm(x)    
plot(x[order(x)],fx[order(x)],yaxt="n",xaxt="n",xlab="",ylab="",type="l",lty=2,col="grey")  
axis(4) 
mtext(side=4,"Density",line=2)  
legend("topleft",c("y", "x density"),   
col=c("black","grey"), lty=1:2, lwd=c(1,2), bty="n")


PDF via method change of variables:    
g1.c<-(-b+sqrt((b^2)-4*a*(c-y)))/(2*a)  
g2.c<-(-b-sqrt((b^2)-4*a*(c-y)))/(2*a)  
g1.prime.c<-abs(1/sqrt((b^2)-4*a*(c-y)))    
fy<-dnorm(g1.c)*abs(g1.prime.c)+    
dnorm(g2.c)*abs(g1.prime.c) 
min(y)  
d<-c+(-(b^2)/(4*a)) 
plot(y,fy,type="l",lwd=2,ylab="density of y",xlab="y", ylim=c(0,0.015), 
main="y=80+30x+6x^2")   
lines(c(44.4,44.4),c(-1,0.01),lty=2)    
lines(c(d,d),c(-1,max(fy)),lty=2,col="red") 
legend("topright", c("d=42.5","d+e=44.4"),lty=2,col=c("red","black")) 


See how it shoots up?? Results are much worse if you choose larger values for the 'a' parameter.
PDF via the CDF method:
d<-c+(-(b^2)/(4*a)) 
first<- 1/(2*sqrt(a)*sqrt(y-d)) 
in_a1<-sqrt(y-d)/sqrt(a)    
in_a2<--sqrt(y-d)/sqrt(a)   
in_b<-b/(2*a)   
A<-in_a1-in_b   
B<-in_a2-in_b   
d   
min(y)  
fy_cdf<-first*(dnorm(A)+dnorm(B))   
plot(y,fy_cdf,type="l",lwd=2,ylab="density of y",xlab="y", ylim=c(0,0.015), 
main="y=80+30x+6x^2")   
lines(c(44.4,44.4),c(-1,0.01),lty=2)    
lines(c(d,d),c(-1,max(fy)),lty=2,col="red") 
legend("topright", c("d=42.5","d+e=44.4"),lty=2,col=c("red","black"))

Note that the results are the same whatever methods is used to derive the pdf: 
# library("miscTools")  
# compPlot(fy_cdf,fy)   
 diff<-fy_cdf-fy   
 summary(abs(diff)) # these are minor rounding errors, I have no issue with that.

An expediment would be to compute either the experimental density or the histogram with:
hist(y)
plot(density(y))

but I would rather have the exact density as I need to compute several statistics from it, notably the probability mass and partial lower mean of some target values. And the derivative of all these with respect to other variables I will interact with X. And well, it is also just for the sack of doing in R what I am writing in the paper...
 A: I will treat the case $a > 0$.  The other cases can be handled with the same methods.
Let $\Phi$ be the CDF of the standard Normal distribution governing the random variable $X$.  Because the graph of $y = ax^2 + bx + c$ is a parabola, the chance that $Y = aX^2 + bX + c$ is less than or equal to $t$ is the chance that $x$ lies between the two roots of $ax^2 + bx + c - t$ (which are easily computed with a quadratic formula).  Calling them $x_{+}(t)$ and $x_{-}(t)$, this is given by $\Phi(x_{+}(t)) - \Phi(x_{-}(t))$.  Differentiation with respect to $t$ gives the PDF:
$$\phi(t) = \frac{1}{\sqrt{2 \pi  d}} e^{-\frac{d}{2 a^2}} \left(e^{\frac{\left(b-\sqrt{d}\right)^2}{8 a^2}}+e^{\frac{\left(b+\sqrt{d}\right)^2}{8 a^2}}\right) e^{-\frac{c-t}{a}}$$
where $d=b^2 - 4a(c-t)$ is the discriminant of the quadratic form and $t$ is not less than the minimum of $y$, equal to $c - b^2/4a.$  Note the factor of $\sqrt{d}$ in the denominator: as $d$ approaches $0$ from above (that is, $c-t$ approaches $b^2/4a$ from higher values of $t$), the exponentials stay finite but the entire PDF has to diverge.  Necessarily this PDF will have a vertical asymptote at $t = c - b^2/4a$. 
Plots of a typical CDF and PDF (for $a=1, b=-2, c=0$} are:

Because these plots coincide with the (obviously correct) plot of a $\chi^2(1)$ distribution when $b=0$--the two plots just get translated one unit to the right--they likely are correct.  Let's look at their counterparts for the example $a=6, b=30, c=80$ in the question:

The PDF looks just like that shown in the question, suggesting the R code actually is correct.
The reason the PDF shoots up at the left is made clear by the CDF: the quadratic behavior of the parabola near its vertex at $-30/12$, where it has a height of $85/2$, focuses a great deal of probability just to the right of $t=85/2$.  The infinite density reflects the verticality of the CDF at that location.  (The verticality is clear in the first plot; in the second it is apparent only when zoomed in closely near its left endpoint.)  Values of $Y$ smaller than $85/2$ cannot be achieved at all, whence there is no probability associated with $t\lt 85/2$.
