# How can I plot all solutions for a system of equations across a given domain for a selected variable?

I have a coupled chemical equilibria system consisting of 12 equations in 12 unknowns. I have managed use the Solve function to write a command that will find a particular solution at any fixed value of pH. I need to plot the value of each unknown across a given pH range. How can I do this?

• Without seeing the code, it's hard to say (you can post the code here, and if it's large, try to condense it to just the most important part). – Jason B. Sep 8 '16 at 21:14

You sure have to wait sometimes to get a chemistry question around here lol. Perhaps this helps you in your PChem studies.

In the absence of your actual equations, let's take a similar example to the one you describe. Perhaps we know the initial concentration of phosphoric acid added , and we want to know the concentrations of all the various species in solution.

The chemical reactions involved are

That's four equations and six unknowns, but you can get two more equations by conservation of phosphate species, and conservation of protons (assuming an initial concentration of $10^{-7}$. You can feed all that to Solve, and if you also give conditions that the concentrations be positive then you get back only one solution. You've done that, but how do you plot it?

When you use Solve, it returns a list of replacement rules, like this

Solve[{x + y == 2, x - y == 2}, {x, y}]
(* {{x -> 2, y -> 0}} *)


To get the values, you need to use ReplaceAll, or it's infix form /.

{x, y} /. Solve[{x + y == 2, x - y == 2}, {x, y}]
(* {{2, 0}} *)


So here is the phosphoric acid concentration calculator:

phosphoricAcidConcentrations[
initConc_] :=
{h, oh, h3PO4, h2PO4, hPO4, pO4} /. Quiet@Solve[
{(h h2PO4)/h3PO4 == kA1,
(h hPO4)/h2PO4 == kA2,
(h pO4)/hPO4 == kA3,
h oh == kW,
h3PO4 + h2PO4 + hPO4 + pO4 == initConc,
h == 3 pO4 + 2 hPO4 + h2PO4 + 10^-7,
h > 0,
h2PO4 > 0,
h3PO4 > 0,
pO4 > 0},
{h, oh, h3PO4, h2PO4, hPO4, pO4}]


(as a PChem student, you'll probably want to include activity coefficients, but I'll leave them out for brevity). Let's test it by taking getting the pH for a 0.1M solution

phosphoricAcidConcentrations[0.1]
-Log10[%]
(* {{0.0235359, 4.24882*10^-13, 0.0764642, 0.0235357,
6.20863*10^-8, 1.18923*10^-18}} *)
(* {{1.62827, 12.3717, 1.11654, 1.62827, 7.207, 17.9247}} *)


You can compare this to what you get from Alpha

WolframAlpha["pH 0.1M phosphoric acid", "Result"]
(* 1.63 *)


Looks good,* so how to plot it? I would collect the results in a list, and use ListLinePlot,

result = Transpose@Table[
First@phosphoricAcidConcentrations[10^(-x)], {x, 0, 10, .1}];

ListLinePlot[-Log10[result],
Frame -> True, FrameLabel -> {"-log(initial conc.)"},
PlotLegends -> {"-log[H]", "-log[OH]",
"-log[\!$$\*SubscriptBox[\(H$$, $$3$$]\)\!$$\*SubscriptBox[\(PO$$, \
$$4$$]\)]",
"-log[\!$$\*SubscriptBox[\(H$$, \
$$2$$]\)\!$$\*SubsuperscriptBox[\(PO$$, $$4$$, $$-$$]\)]",
"-log[\!$$\*SubsuperscriptBox[\(HPO$$, $$4$$, $$\(2$$$$-$$\)]\)]",
"-log[\!$$\*SubsuperscriptBox[\(PO$$, $$4$$, $$\(3$$$$-$$\)]\)]"}]


Excuse the subscript-box gobbledegook.

*: But don't try it with very low concentrations, Alpha doesn't treat this as a polyprotic weak acid I think.