The functions Efunction and dEfunction compute the expectation value of the energy and its derivative. They use the the quasi-Newton method of Broyden, Fletcher, Goldfarb, and Shanno (BFGS) It uses the first derivatives only. The BFGS algorithm has proven good performance even for non-smooth optimizations. These functions need to be changed when you want to your own derivatives.

//  this function defines the expectation value of the local energy
double Efunction(Vector  &x)
{
  double value = x(0)*x(0)*0.5+1.0/(8*x(0)*x(0));
  return value;
} // end of function to evaluate

//  this function defines the derivative of the energy 
void dEfunction(Vector &x, Vector &g)
{
  g(0) = x(0)-1.0/(4*x(0)*x(0)*x(0));
} // end of function to evaluate

You need to change these functions in order to compute the local energy for your system. I used 1000 cycles per call to get a new value of \( \langle E_L[\alpha]\rangle \). When I compute the local energy I also compute its derivative. After roughly 10-20 iterations I got a converged result in terms of \( \alpha \).