diff --git a/include/cantera/base/ct_defs.h b/include/cantera/base/ct_defs.h
index 83717fc45..f8014d778 100644
--- a/include/cantera/base/ct_defs.h
+++ b/include/cantera/base/ct_defs.h
@@ -164,7 +164,13 @@ const doublereal MaxExp = 690.775527898;
 //! Fairly random number to be used to initialize variables against
 //! to see if they are subsequently defined.
 const doublereal Undef = -999.1234;
+
 //! Small number to compare differences of mole fractions against.
+/*!
+ *   This number is used for the interconversion of mole fraction and mass fraction quantities
+ *   when the molecuar weight of a species is zero. It's also used for the matrix inversion
+ *   of transport properties when mole fractions must be positive.
+ */
 const doublereal Tiny = 1.e-20;
 
 //! inline function to return the max value of two doubles.
diff --git a/include/cantera/thermo/Phase.h b/include/cantera/thermo/Phase.h
index d632a4bf5..d4cbc50c0 100644
--- a/include/cantera/thermo/Phase.h
+++ b/include/cantera/thermo/Phase.h
@@ -9,7 +9,7 @@
 
 #include "cantera/base/vec_functions.h"
 #include "cantera/base/ctml.h"
-#include "Elements.h"
+#include "cantera/thermo/Elements.h"
 
 namespace Cantera
 {
diff --git a/src/transport/Tortuosity.h b/include/cantera/transport/Tortuosity.h
similarity index 100%
rename from src/transport/Tortuosity.h
rename to include/cantera/transport/Tortuosity.h
diff --git a/src/numerics/BandMatrix.cpp b/src/numerics/BandMatrix.cpp
index bf41b1699..096d4872c 100644
--- a/src/numerics/BandMatrix.cpp
+++ b/src/numerics/BandMatrix.cpp
@@ -160,7 +160,9 @@ doublereal BandMatrix::value(size_t i, size_t j) const
 //====================================================================================================================
 size_t BandMatrix::index(size_t i, size_t j) const
 {
-    size_t rw = m_kl + m_ku + i - j;
+    int jj = j;
+    int ii = i;
+    size_t rw = (int) m_kl + (int) m_ku + (int) ii - jj;
     return (2*m_kl + m_ku + 1)*j + rw;
 }
 //====================================================================================================================
@@ -218,12 +220,14 @@ vector_int&   BandMatrix::ipiv()
  */
 void BandMatrix::mult(const doublereal* b, doublereal* prod) const
 {
-    size_t nr = nRows();
+    int kl = m_kl;
+    int ku = m_ku;
+    int nr = nRows();
     doublereal sum = 0.0;
-    for (size_t m = 0; m < nr; m++) {
+    for (int m = 0; m < nr; m++) {
         sum = 0.0;
-        for (size_t j = m - m_kl; j <= m + m_ku; j++) {
-            if (j < m_n) {
+        for (int j = m - kl; j <= m + ku; j++) {
+            if (j >= 0 && j < (int) m_n) {
                 sum += _value(m,j) * b[j];
             }
         }
@@ -236,13 +240,16 @@ void BandMatrix::mult(const doublereal* b, doublereal* prod) const
  */
 void BandMatrix::leftMult(const doublereal* const b, doublereal* const prod) const
 {
-    size_t nc = nColumns();
+    int kl = m_kl;
+    int ku = m_ku;
+    int nc = nColumns();
     doublereal sum = 0.0;
-    for (size_t n = 0; n < nc; n++) {
+    for (int n = 0; n < nc; n++) {
         sum = 0.0;
-        for (size_t i = n - m_ku; i <= n + m_kl; i++) {
-            if (i < m_n) {
-                sum += _value(i,n) * b[i];
+        for (int i = n - ku; i <= n + kl; i++) {
+            if (i >= 0 && i < (int) m_n) {
+                size_t ii = i;
+                sum += _value(ii,n) * b[ii];
             }
         }
         prod[n] = sum;
@@ -421,12 +428,14 @@ int BandMatrix::factorAlgorithm() const
 // Returns the one norm of the matrix
 doublereal BandMatrix::oneNorm() const
 {
+    int ku = m_ku;
+    int kl = m_kl;
     doublereal value = 0.0;
-    for (size_t j = 0; j < m_n; j++) {
+    for (int j = 0; j < (int) m_n; j++) {
         doublereal sum = 0.0;
         doublereal* colP =  m_colPtrs[j];
-        for (size_t i = j - m_ku; i <= j + m_kl; i++) {
-            sum += fabs(colP[m_kl + m_ku + i - j]);
+        for (int i = j - ku; i <= j + kl; i++) {
+            sum += fabs(colP[kl + ku + i - j]);
         }
         if (sum > value) {
             value = sum;
@@ -440,10 +449,10 @@ size_t BandMatrix::checkRows(doublereal& valueSmall) const
     valueSmall = 1.0E300;
     size_t iSmall = npos;
     double vv;
-    for (size_t i = 0; i < m_n; i++) {
+    for (int i = 0; i < (int) m_n; i++) {
         double valueS = 0.0;
-        for (size_t j = i - m_kl; j <= i + m_ku; j++) {
-            if (j < m_n) {
+        for (int j = i - (int) m_kl; j <= i + (int) m_ku; j++) {
+            if (j >= 0 && j < (int) m_n) {
                 vv = fabs(value(i,j));
                 if (vv > valueS) {
                     valueS = vv;
@@ -466,10 +475,10 @@ size_t BandMatrix::checkColumns(doublereal& valueSmall) const
     valueSmall = 1.0E300;
     size_t jSmall = npos;
     double vv;
-    for (size_t j = 0; j < m_n; j++) {
+    for (int j = 0; j < (int) m_n; j++) {
         double valueS = 0.0;
-        for (size_t i = j - m_ku; i <= j + m_kl; i++) {
-            if (i < m_n) {
+        for (int i = j - (int) m_ku; i <= j + (int) m_kl; i++) {
+            if (i >= 0 && i < (int) m_n) {
                 vv = fabs(value(i,j));
                 if (vv > valueS) {
                     valueS = vv;
diff --git a/src/numerics/DAE_solvers.cpp b/src/numerics/DAE_solvers.cpp
new file mode 100644
index 000000000..a4b9cc7e1
--- /dev/null
+++ b/src/numerics/DAE_solvers.cpp
@@ -0,0 +1,45 @@
+/**
+ *  @file DAE_solvers.cpp
+ *       Factory routine for picking the DAE solver package
+ */
+/*
+ * $Revision: 725 $
+ * $Date: 2011-05-16 18:45:08 -0600 (Mon, 16 May 2011) $
+ */
+/*
+ * Copyright 2004 Sandia Corporation. Under the terms of Contract
+ * DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government
+ * retains certain rights in this software.
+ * See file License.txt for licensing information.
+ */
+
+
+#include "cantera/base/ct_defs.h"
+#include "cantera/numerics/DAE_Solver.h"
+#include "cantera/numerics/IDA_Solver.h"
+
+// DAE_DEVEL is turned off at the current time
+#define DAE_DEVEL
+#ifdef DAE_DEVEL
+
+
+namespace Cantera {
+
+  DAE_Solver* newDAE_Solver(std::string itype, ResidJacEval& f) {
+        if (itype == "IDA") {
+#ifdef HAS_SUNDIALS
+            return new IDA_Solver(f);
+#else
+            throw CanteraError("newDAE_Solver","IDA solver requires sundials"
+			       " package, but Cantera was not built with sundials.");
+#endif
+        }
+        else {
+            throw CanteraError("newDAE_Solver",
+			       "unknown DAE solver: "+itype);
+        }
+    }
+}
+
+#
+#endif
diff --git a/src/numerics/NonlinearSolver.cpp b/src/numerics/NonlinearSolver.cpp
index 8925523a2..72cc8d594 100644
--- a/src/numerics/NonlinearSolver.cpp
+++ b/src/numerics/NonlinearSolver.cpp
@@ -754,7 +754,6 @@ void NonlinearSolver::scaleMatrix(GeneralMatrix& jac, doublereal* const y_comm,
                                   doublereal time_curr, int num_newt_its)
 {
     size_t irow, jcol;
-    size_t ku, kl;
     size_t ivec[2];
     jac.nRowsAndStruct(ivec);
     double* colP_j;
@@ -783,12 +782,12 @@ void NonlinearSolver::scaleMatrix(GeneralMatrix& jac, doublereal* const y_comm,
                     }
                 }
             } else if (jac.matrixType_ == 1) {
-                kl = ivec[0];
-                ku = ivec[1];
-                for (jcol = 0; jcol < neq_; jcol++) {
+                int kl = ivec[0];
+                int ku = ivec[1];
+                for (int jcol = 0; jcol < (int) neq_; jcol++) {
                     colP_j = (doublereal*) jac.ptrColumn(jcol);
-                    for (irow = jcol - ku; irow <= jcol + kl; irow++) {
-                        if (irow < neq_) {
+                    for (int irow = jcol - ku; irow <= jcol + kl; irow++) {
+                        if (irow >= 0 && irow < (int) neq_) {
                             colP_j[kl + ku + irow - jcol] *= m_colScales[jcol];
                         }
                     }
@@ -828,12 +827,12 @@ void NonlinearSolver::scaleMatrix(GeneralMatrix& jac, doublereal* const y_comm,
                 }
             }
         } else if (jac.matrixType_ == 1) {
-            kl = ivec[0];
-            ku = ivec[1];
-            for (jcol = 0; jcol < neq_; jcol++) {
+            int kl = ivec[0];
+            int ku = ivec[1];
+            for (int jcol = 0; jcol < (int) neq_; jcol++) {
                 colP_j = (doublereal*) jac.ptrColumn(jcol);
-                for (irow = jcol - ku; irow <= jcol + kl; irow++) {
-                    if (irow < neq_) {
+                for (int irow = jcol - ku; irow <= jcol + kl; irow++) {
+                    if (irow >= 0 && irow < (int) neq_) {
                         double vv = fabs(colP_j[kl + ku + irow - jcol]);
                         if (m_rowScaling) {
                             m_rowScales[irow] += vv;
@@ -871,12 +870,12 @@ void NonlinearSolver::scaleMatrix(GeneralMatrix& jac, doublereal* const y_comm,
                     }
                 }
             } else if (jac.matrixType_ == 1) {
-                kl = ivec[0];
-                ku = ivec[1];
-                for (jcol = 0; jcol < neq_; jcol++) {
+                int kl = ivec[0];
+                int ku = ivec[1];
+                for (int jcol = 0; jcol < (int) neq_; jcol++) {
                     colP_j = (doublereal*) jac.ptrColumn(jcol);
-                    for (irow = jcol - ku; irow <= jcol + kl; irow++) {
-                        if (irow < neq_) {
+                    for (int irow = jcol - ku; irow <= jcol + kl; irow++) {
+                        if (irow >= 0 && irow < (int) neq_) {
                             colP_j[kl + ku + irow - jcol] *= m_rowScales[irow];
                         }
                     }
@@ -3881,14 +3880,24 @@ int NonlinearSolver::beuler_jac(GeneralMatrix& J, doublereal* const f,
 
                 doublereal diff;
 
-
-
+                int ileft =  (int) j - (int) ku;
+                int iright=   j + kl;
+                for (int i = ileft; i <= iright; i++) {
+                    if (i >= 0 &&  i < (int) neq_) {
+                        size_t ii = i;
+                        size_t index = (int) kl + (int) ku + i - (int) j; 
+                        diff = subtractRD(m_wksp[ii], f[ii]);
+                        col_j[index] = diff / dy;
+                    }
+                }
+/*
                 for (size_t i = j - ku; i <= j + kl; i++) {
                     if (i < neq_) {
                         diff = subtractRD(m_wksp[i], f[i]);
                         col_j[kl + ku + i - j] = diff / dy;
                     }
                 }
+*/
                 y[j] = ysave;
                 if (solnType_ != NSOLN_TYPE_STEADY_STATE) {
                     ydot[j] = ydotsave;
diff --git a/src/thermo/Phase.cpp b/src/thermo/Phase.cpp
index edbdad34a..cabf21b33 100644
--- a/src/thermo/Phase.cpp
+++ b/src/thermo/Phase.cpp
@@ -318,6 +318,9 @@ void Phase::setMoleFractions(const doublereal* const x)
 {
     // Use m_y as a temporary work vector for the non-negative mole fractions
     doublereal norm = 0.0;
+    /*
+     * sum is calculated below as the unnormalized molecular weight
+     */
     doublereal sum = 0;
     for (size_t k = 0; k < m_kk; k++) {
         double xk = std::max(x[k], 0.0); // Ignore negative mole fractions
@@ -325,13 +328,19 @@ void Phase::setMoleFractions(const doublereal* const x)
         norm += xk;
         sum += m_molwts[k] * xk;
     }
-
-    transform(m_y.begin(), m_y.end(), m_ym.begin(),
-              timesConstant<double>(1.0/sum));
-
-    // Now set m_y to the mass fractions
-    transform(m_ym.begin(), m_ym.begin() + m_kk, m_molwts.begin(),
-              m_y.begin(), multiplies<double>());
+    /*
+     * Set m_ym_ to the normalized mole fractions divided by the normalized mean molecular weight:
+     *         m_ym_k = X_k / (sum_k X_k M_k)
+     */
+    transform(m_y.begin(), m_y.end(), m_ym.begin(), timesConstant<double>(1.0/sum));
+    /*
+     * Now set m_y to the normalized mass fractions
+     *          m_y =  X_k M_k / (sum_k X_k M_k)
+     */
+    transform(m_ym.begin(), m_ym.begin() + m_kk, m_molwts.begin(), m_y.begin(), multiplies<double>());
+    /*
+     * Calculate the normalized molecular weight
+     */
     m_mmw = sum/norm;
 
     // Call a routine to determine whether state has changed.
@@ -961,7 +970,7 @@ void Phase::init(const vector_fp& mw)
         }
 
         // Some surface phases may define species representing empty sites
-        // hat have zero molecular weight. Give them a very small molecular
+        // that have zero molecular weight. Give them a very small molecular
         // weight to avoid dividing by zero.
         if (m_molwts[k] < Tiny) {
             m_molwts[k] = Tiny;