
d_psi_d_delta = -2.0 * y[i].Ci * dm1 * psi;
d2_psi_d_delta2 = (2.0 * y[i].Ci * d2 - 1.0) * 2.0
* y[i].Ci * psi;
d_Delta_d_delta = dm1 * (y[i].Ai * theta * (2.0 /
y[i].betai)
* pow(d2, (0.5/y[i].betai) - 1)
+ 2.0 * y[i].Bi * y[i].ai * pow(d2,y[i].
ai - 1.0));
d = (y[i].Ai / y[i].betai) * pow (d2, 0.5/y[i].
betai - 1.0);
d2_Delta_d_delta2 = (d_Delta_d_delta / dm1)
+ d2 * (4.0 * y[i].Bi * y[i].ai * (y[i].ai - 1)
* pow (d2, y[i].ai - 2.0)
+ 2.0 * d * d
+ y[i].Ai * theta * (4.0/y[i].betai) * (0.5/
y[i].betai - 1.0)
* pow (d2, 0.5/y[i].betai - 2.0));
Delta_bim1 = pow ...