1. No guidance is provided as to the appropriate choice of P and Q by this method
  2. A different matrix must be constructed and solved for each value of P and Q under consideration
  3. There is no way to restrict the range of values and/or signs of the calculated coefficients
The third drawback is the most severe, as divergences  in Eq. \ref{eq:Pade_zero_index} will occur whenever the denominator term  \(\sum_{q=0}^Q\gamma_{i,q}B^q=0\) (causing \(c_i\left(B\right)\ \rightarrow\infty\) ).   If some of the \(\gamma_i\) coefficients provided by the matrix method are negative, the resulting fit may lead to divergences at (real) values of \(B\) within the experimental range of the experiment, spoiling the calculation of \(R\left(T,B\right)\).
When we tried this matrix method for P = Q = 2 for the CX1010 thermometer data shown above, it returned values that resulted in divergences in all but the first two coefficients. We did not investigate the method further (for example, for P = 3 and Q = 1). In comparison, note that divergences cannot occur at real values of \(B\) for  the Padé fits to  \(y_i\left(B\right)\) constructed using our original method — in which we modeled our choice of P and Q  and  the initial values on the Padé expansion of Eq. \ref{eq:exponential_decay_from_zero} — because the Padé expansions of \(e^{-z}\)and \(e^{-z}-1\) only have positive terms in the denominator. Given positive initial values for \(\gamma_{i,q}\), we were then able to constrain the fit to only consider positive values through the specification of appropriate boundary_values.  For further details, see the discussion regarding the setting of parameter boundaries in the  Scipy manual entry for curve_fit .