[R] solving nonlinear equations in R

[Demirtas, Hakan]

Hi, I can't seem to get computationally stable estimates for the following
system:

Y=a+bX+cX^2+dX^3, where X~N(0,1). (Y is expressed as a linear combination
of the first three powers of a standard normal variable.) Assuming that
E(Y)=0 and Var(Y)=1, one can obtain the following equations after tedious
algebraic calculations:

1) b^2+6bd+2c^2+15d^2=1
2) 2c(b^2+24bd+105d^2+2)=E(Y^3)
3) 24[bd+c^2(1+b^2+28bd)+d^2(12+48bd+141c^2+225d^2)]=E(Y^4)-3

Obviously, a=-c. Suppose that distributional form of Y is given so we know
E(Y^3) and E(Y^4). In other words, we have access to the third and fourth
raw moments. How do we solve for these four coefficients? I reduced the
number of unknowns/equations to two, and subsequently used a grid
approach. It works well when I am close to the center of the support, but
fails miserably at the tails. Any ideas? Hopefully, there is a nice R
function that does this.

Hakan Demirtas