It's an eiganvalue problem. See the normal roots implementation here.

https://github.com/numpy/numpy/blob/master/numpy/polynomial/polynomial.py

The implicit smoother has a large deviation from the input function at the edges.

Upgrade data structures from json to hdf5.

charts can load json with format {name:data for name in names}

Do this for properties with property_name:data, for concentrations with ion:data.

Use matrix math as in spresso M.full.

Each time point should be represented by a data structure in the same format as the initialization data structure.

Consider both nodes and concentrations.

Rounding error in the derivatives cause the flat function to have small errors in derivative. When these are divided by the maximum of that same derivative, they are the same order as other derivatives.

• Equilibration base.
• Flux base.

pH solver uses cH^n, where n is number of valence species.

perhaps perform pH calculations in log space?