There may be many reasons for "errors".
For example you provide only a small number of decimals for your constants (4 - 5) while sometimes you provide 16 decimals. It might be better to provide at least 16 for all (and more if it turms out that higher precision is used) to reduce rounding errors
This is a bit related to a mathematical fact: generically you can expect discrete solutions if using n=13 equations and n variables (discrete = isolated n-tuples in n-space each giving a solution). But even if there would be only discrete solutions there might be none, especially over the Reals (think pf a parabola).
And if feeding the solutions there will be a (numerical) error. One has to judge whether it can be accepted or not - for example a very steep function with a zero z0 may have a large absoluet error in that while it can be accepted because it is small in a relative sense (and the "best" approximation).
From you last eq13 and your remark "must be positive" it follows that u_C1 is between ionic=0 and ionic=infinity, i.e. roughly between 0.58 and 0.97