diff --git a/include/teqp/derivs.hpp b/include/teqp/derivs.hpp
index 51dd9506df418875efcc98efb9a239dadcbe07c4..6d280701ce6466a5e7a1d1068108eddfcda9d507 100644
--- a/include/teqp/derivs.hpp
+++ b/include/teqp/derivs.hpp
@@ -3,6 +3,7 @@
 #include <optional>
 #include <complex>
 #include <tuple>
+#include <map>
 
 #include "MultiComplex/MultiComplex.hpp"
 
@@ -110,13 +111,27 @@ typename ContainerType::value_type get_B2vir(const Model& model, const TType T,
     return B2;
 }
 
+/*
+ B_n = \frac{1}{(n-2)!} lim_rho\to 0 d^{n-1}alphar/drho^{n-1}|T,z
+*/
+
 template <typename Model, typename TType, typename ContainerType>
 auto get_Bnvir(const Model& model, int Nderiv, const TType T, const ContainerType& molefrac) {
-    // B_n = lim_rho\to 0 d^{n-1}alphar/drho^{n-1}|T,z
+    
     using namespace mcx;
     using fcn_t = std::function<MultiComplex<double>(const MultiComplex<double>&)>;
     fcn_t f = [&model, &T, &molefrac](const MultiComplex<double>& rho_) -> MultiComplex<double> { return model.alphar(T, rho_, molefrac); };
-    return diff_mcx1(f, 0.0, Nderiv, true /* and_val */);
+    std::map<int, TType> o;
+    auto dalphardrhon = diff_mcx1(f, 0.0, Nderiv+1, true /* and_val */);
+    for (int n = 2; n < Nderiv+1; ++n) {
+        o[n] = dalphardrhon[n-1];
+        // 0!=1, 1!=1, so only n>3 terms need factorial correction
+        if (n > 3) {
+            auto factorial = [](int N) {return tgamma(N + 1); };
+            o[n] /= factorial(n-2);
+        }
+    }
+    return o;
 }
 
 /***
diff --git a/src/main.cpp b/src/main.cpp
index e7447cb92a94c326100c5439bd58357424541beb..1ef6bbae864596ddead266770810d0ac725bd561 100644
--- a/src/main.cpp
+++ b/src/main.cpp
@@ -42,17 +42,17 @@ void test_vdW() {
     auto Nvir = 8;
     auto Bn = get_Bnvir(vdW, Nvir, T, molefrac);
     // Exact solutions for virial coefficients for van der Waals 
-    auto get_vdW_exacts = [a,b,R,T](int additional){ 
-        std::vector<double> o = {0, b - a / (R * T)};
-        for (auto i = 0; i < additional; ++i) {
-            auto factorial = [](int N){return tgamma(N + 1);};
-            o.push_back(factorial(i+1)*pow(b,i+2));
+    auto get_vdW_exacts = [a,b,R,T](int Nmax){
+        std::map<int, double> o = {{2, b - a / (R * T)}};
+        for (auto i = 3; i <= Nmax; ++i) {
+            o[i] = pow(b, i-1);
         }
-        return std::valarray<double>(&(o[0]), o.size());
+        return o;
     };
-    std::valarray<double> Bnexact = get_vdW_exacts(Nvir-1);
-    std::valarray<double> errrr = Bnexact - std::valarray<double>(&(Bn[0]), Bn.size());
-
+    auto Bnexact = get_vdW_exacts(Nvir);
+    for (auto i = 2; i <= Nvir; ++i){
+        std::cout << std::scientific << i << ", " << Bnexact[i] << ", " << Bn[i] << ", " << std::abs(Bnexact[i]-Bn[i]) << std::endl;
+    }
     int rr = 0;
 
 }