# Making a Complicated Useable Function using NSolve

I have been working on this problem for awhile and cannot seem to get a good answer. Essentially, I want to plot existence regions with RegionPlot using equations for parameters found using NSolve. Here is the code:

gB = .15;
gR = .15;
MR = .5;
MB = .5;

{mGB, mBB, mRR, mBO, mRO, mBR, mGBO} =
NSolve[mRR == MR - ((10 - nr)*(1 - gB)/nr)*mBR && mRO == 1 - mRR &&
mGB == ((MB - (1 - gB)*mBB)/(gB)) && mRR*mBO == mRO*mBR &&
mBO*mGB == mBB*mGBO && mBO == 1 - mBB - mBR && mGBO == 1 - mGB &&
mGB >= 0 && mBB >= 0 && mRR >= 0 && mBO >= 0 && mRO >= 0 &&
mBR >= 0 && mGBO >= 0 && nr >= 0 && nr <= 10.000001, {mGB, mBB,
mRR, mBO, mRO, mBR, mGBO}]


So I am solving for all the mII as a function of nr. I would then like to create a function using this output for each mII. I would then like to be able to call these functions using RegionPlot to find the existence region for other inequalities that depend on these values while varying the value of nr. I have tried ReplaceAll with no success. I am fairly new to Mathematica, so I suspect there is something going on that I do not understand. Hopefully this is clear!

You can even solve exactly, if you rationalize.

(li = List @@
Rationalize[
mRR == MR - ((10 - nr)*(1 - gB)/nr)*mBR && mRO == 1 - mRR &&
mGB == ((MB - (1 - gB)*mBB)/(gB)) && mRR*mBO == mRO*mBR &&
mBO*mGB == mBB*mGBO && mBO == 1 - mBB - mBR && mGBO == 1 - mGB &&
mGB >= 0 && mBB >= 0 && mRR >= 0 && mBO >= 0 && mRO >= 0 &&
mBR >= 0 && mGBO >= 0 && 0 <= nr <= 10.000001, 0])

sol = FullSimplify[
Solve[li, {mGB, mBB, mRR, mBO, mRO, mBR, mGBO}, Method -> Reduce],
0 < nr < 10000001/1000000]

(*    {{mGB -> (55 + 11 nr - Sqrt[7225 + 4 nr (340 + nr)])/(-60 + 9 nr),
mBB -> (-255 + 19 nr + Sqrt[7225 + 4 nr (340 + nr)])/(-340 + 51 nr),
mRR -> (nr (5 + 28 nr - Sqrt[7225 + 4 nr (340 + nr)]) +
10 (-85 + Sqrt[7225 + 4 nr (340 + nr)]))/(2 nr (70 + 13 nr)),
mBO -> (11050 - 270 Sqrt[7225 + 4 nr (340 + nr)] +
nr (1585 - 34 nr + 17 Sqrt[7225 + 4 nr (340 + nr)]))/(
17 (-20 + 3 nr) (70 + 13 nr)),
mRO -> (-10 (-85 + Sqrt[7225 + 4 nr (340 + nr)]) +
nr (135 - 2 nr + Sqrt[7225 + 4 nr (340 + nr)]))/(
2 nr (70 + 13 nr)),
mBR -> -((10 (-85 - 15 nr + Sqrt[7225 + 4 nr (340 + nr)]))/(
17 (70 + 13 nr))),
mGBO -> (-115 - 2 nr + Sqrt[7225 + 4 nr (340 + nr)])/(-60 + 9 nr)}}    *)

Needs["PlotLegends"]

Plot[Evaluate[{mGB, mBB, mRR, mBO, mRO, mBR, mGBO} /. First@sol], {nr,
0, 10000001/1000000},
PlotStyle -> Table[Hue[.8 (i - 1)/6], {i, 1, 7}],
PlotLegend -> {mGB, mBB, mRR, mBO, mRO, mBR, mGBO},
LegendPosition -> {1.1, -0.4}, ImageSize -> 500] (min = Minimize[{#, 0 < nr < 10000001/1000000}, nr] & /@
Evaluate[{mGB, mBB, mRR, mBO, mRO, mBR, mGBO} /. First@sol] //
Quiet) // N

(*    {{0.5, {nr -> 0.}}, {0.474635, {nr -> 10.}}, {0., {nr ->
0.}}, {0.262682, {nr -> 10.}}, {0.5, {nr -> 10.}}, {0., {nr ->
0.}}, {0.356267, {nr -> 10.}}}    *)

(max = Maximize[{#, 0 < nr < 10000001/1000000}, nr] & /@
Evaluate[{mGB, mBB, mRR, mBO, mRO, mBR, mGBO} /. First@sol] //
Quiet) // N

(*    {{0.643733, {nr -> 10.}}, {0.5, {nr -> 0.}}, {0.5, {nr ->
10.}}, {0.5, {nr -> 0.}}, {1., {nr -> 0.}}, {0.262682, {nr ->
10.}}, {0.5, {nr -> 0.}}}    *)
`