Here is a corrected version of your code with a few notes: #include <iostream> #include <complex> #include <boost/foreach.hpp> using namespace std; // use rational numbers #include <boost/rational.hpp> // use math constants #include <boost/math/constants/constants.hpp> // use math factorial #include <boost/math/special_functions/factorials.hpp> //using namespace boost::math // use units #include <boost/typeof/std/complex.hpp> // MCS - need to include this header for make_dimension_list #include <boost/mpl/list.hpp> #include <boost/units/cmath.hpp> // MCS - need to include this header #include <boost/units/pow.hpp> #include <boost/units/io.hpp> #include <boost/units/systems/si.hpp> #include <boost/units/systems/si/prefixes.hpp> #include <boost/units/systems/si/codata/electron_constants.hpp> #include <boost/units/systems/si/codata/physico-chemical_constants.hpp> #include <boost/units/systems/si/codata/universal_constants.hpp> using namespace boost::units; using namespace boost::units::si; // MCS - need codata namespace to find hbar using namespace boost::units::si::constants::codata; //use units end vector<boost::rational<signed long> > hermite_polynomial_coefficients(int order, boost::rational<signed long> lambda_per_alpha) { vector<boost::rational<signed long> > c_even; // can't use static, because lambda_per_alpha changes between calls vector<boost::rational<signed long> > c_odd; if(c_even.size() == 0) { c_even.push_back(1); c_even.push_back(0); } if(c_odd.size() == 0) { c_odd.push_back(0); c_odd.push_back(1); } vector<boost::rational<signed long> >& a( (order % 2 == 0) ?c_even:c_odd ); if(order+1 < a.size()) return a; for(int i = a.size() ; i <= order ; ++i) { boost::rational<signed long> next = a[i-2]*(2*(i-2)+1-lambda_per_alpha)/(i*(i-1)); a.push_back(next); } // BOOST_FOREACH(boost::rational<signed long> v,a) // std::cout << v << " "; // std::cout << "\n"; return a; } vector<boost::rational<signed long> > hermite_polynomial_scaled(int order, boost::rational<signed long> lambda_per_alpha) { vector<boost::rational<signed long> > a(hermite_polynomial_coefficients(order,lambda_per_alpha)); vector<boost::rational<signed long> > b; boost::rational<signed long> factor(std::pow(2.,order)/a[order]); for(int i = 0 ; i<=order ; ++i) b.push_back(a[i]*factor); return b; } // MCS - derived_dimension only works with integer powers; use make_dimension_list here instead typedef make_dimension_list<boost::mpl::list<dim<length_base_dimension,static_rational<-1,2> > > >::type quantum_wavefunction_1D_dimension; typedef unit<quantum_wavefunction_1D_dimension,si::system> quantum_wavefunction_1D; // this function returns 1/sqrt(meter) quantity<quantum_wavefunction_1D,complex<double> > quantum_oscillator_wavefunction( int n // n - order of wavefunction , const quantity<length> x // position , const quantity<mass> mass // mass - mass of particle , const quantity<frequency> omega // oscillator frequency ) { // MCS - need to define unit typedef in addition to dimension typedef derived_dimension<length_base_dimension,-2>::type reciprocal_area_dimension; typedef unit<reciprocal_area_dimension,si::system> reciprocal_area; quantity<energy_time> h_(hbar); quantity<energy> E(h_*omega*(n+0.5)); // energy of oscillator quantity<reciprocal_area> lambda(2.0*mass*E/(h_*h_)), alpha(mass*omega/h_); quantity<wavenumber> alpha_root(root<2>(alpha)); double lambda_per_alpha(2*n+1); const double pi(boost::math::constants::pi<double>()), fac(boost::math::factorial<double>(n)); // N_n should be 1/sqrt(meter) because it has a dimension of quantum_wavefunction_1D const quantity<quantum_wavefunction_1D,complex<double> > N_n(root<2>(root<2>(alpha/pi)/(std::pow(2.,n)*fac))); const complex<double> i(0,1); quantity<quantum_wavefunction_1D,complex<double> > result(0/root<2>(meter)); vector<boost::rational<signed long> > a(hermite_polynomial_scaled(n,lambda_per_alpha)); for(int v=0 ; v<=n ; v++) { // result += N_n*a[v]*pow<v>(alpha_root*x)*exp(-0.5*alpha*pow<2>(x)); // MCS - a[v] is boost::rational -> need to convert to double via rational_cast // MCS - pow<v> does not work because v must be known at compile time const double a_v(boost::rational_cast<double>(a[v])), tmp1(a_v*pow(alpha_root*x,v)); const complex<double> tmp2(exp(-0.5*alpha*pow<2>(x))), tmp3(tmp1*tmp2); result += tmp3*N_n; } return result; } int main() { // MCS - change -1 -> +1 here otherwise loop doesn't run quantity<length> s(15*nano*meter), ds(0.01*nano*meter); std::cout << s << "\t" << ds << std::endl; for (quantity<length> x=-s ; x<=s ; x+=ds) { const quantity<temperature> T(273.15*kelvin); const quantity<frequency> omega(k_B*T/hbar.value()); std::cout << x << " " << quantum_oscillator_wavefunction(25,x,m_e,omega) << "\n"; } for (quantity<length> x=-s ; x<=s ; x+=ds) { const quantity<temperature> T(273.15*kelvin); const quantity<frequency> omega(k_B*T/hbar.value()); const quantity<quantum_wavefunction_1D,complex<double> > psi(quantum_oscillator_wavefunction(25,x,m_e,omega)); std::cout << x.value() << " " << psi.value().real() << "\t" << psi.value().imag() << "\n"; } return 0; } and the n=25 wavefunction: Cheers, Matthias