# Numerically solve system of equations with nested variables

I'm trying to solve the following system of equations. The first two equations depend on the variable TeV

RxA[TeV_] :=   Exp[E1/TeV]/ngas*(2*Pi*(me*mp/mp)*qe*TeV/hp^2)^(3/2);(*kT in eV*)
RxB[TeV_] :=  Exp[E2/TeV]/ngas*(2*Pi*(mp^2/(2*mp))*qe*TeV/hp^2)^(3/2);(*kT in eV*)
eqpart = NSolve[{2*y^2/(x - y) == RxA,4*(x - y)^2/(1 - x) == RxB}, {x, y}];
Plot[Evaluate[{x, y} /. eqpart], {TeV, 0.1, 10}, PlotRange -> All]


I don't understand why there is no output on the plot. Is the syntax ok?

Clear["Global*"]


You need to provide values for all of the parameters. I have arbitrarily set them all to 1

me = mp = qe = hp = E1 = E2 = ngas = 1;

RxA[TeV_] := Exp[E1/TeV]/ngas*(2*Pi*(me*mp/mp)*qe*TeV/hp^2)^(3/2);(*kT in eV*)
RxB[TeV_] :=
Exp[E2/TeV]/ngas*(2*Pi*(mp^2/(2*mp))*qe*TeV/hp^2)^(3/2);(*kT in eV*)


You must include arguments for RxA and RxB in the equations, and since eqpart uses a numeric technique it should only be called when TeV has a numeric value.

eqpart[TeV_?NumericQ] :=
NSolve[{2*y^2/(x - y) == RxA[TeV], 4*(x - y)^2/(1 - x) == RxB[TeV]}, {x, y}];

Show[
Plot[#[[1]] /. eqpart[TeV], {TeV, 0.1, 10},
PlotPoints -> 50,
PlotStyle -> ColorData[97][#[[2]]]] & /@
{{x, 1}, {y, 2}},
Frame -> True,
Axes -> False,
FrameLabel -> (Style[#, 14, Bold] & /@ {"TeV", "x, y"}),
PlotRange -> All]


Or, if you want a ParametricPlot

ParametricPlot[{x, y} /. eqpart[TeV],
{TeV, 0.1, 10}, AspectRatio -> 1]


EDIT: To use Solve rather than NSolve (and again setting all parameters to 1)

eqns = {2*y^2/(x - y) == RxA[TeV],
4*(x - y)^2/(1 - x) == RxB[TeV]}

(* {(2 y^2)/(x - y) == 2 Sqrt[2] E^(1/TeV) π^(3/2) TeV^(3/2), (4 (x - y)^2)/(
1 - x) == E^(1/TeV) π^(3/2) TeV^(3/2)} *)


By assuming that x != y && x != 1 the equations can be put in a more manageable form.

eqns2 = eqns // Simplify[#, x != y && x != 1] &

(* {Sqrt[2] E^(1/TeV) π^(3/2) TeV^(3/2) (x - y) == y^2,
E^(1/TeV) π^(3/2) TeV^(3/2) (-1 + x) + 4 (x - y)^2 == 0} *)


The solutions are each a ConditionalExpression with the condition that TeV > 0. Simplifying with this assumption eliminates the condition. The solutions are expressed as Root objects.

sol = Solve[eqns2, {x, y}, Reals] //
Simplify[#, TeV > 0] &

(* {{x -> Root[
E^(4/TeV) π^6 TeV^6 -
64 E^(6/TeV) π^9 TeV^9 + (-4 E^(4/TeV) π^6 TeV^6 -
32 Sqrt[2] E^(4/TeV) π^6 TeV^6 +
128 E^(6/TeV) π^9 TeV^9) #1 + (390 E^(4/TeV) π^6 TeV^6 +
96 Sqrt[2] E^(4/TeV) π^6 TeV^6 -
64 E^(6/TeV) π^9 TeV^9) #1^2 + (-772 E^(4/TeV) π^6 TeV^6 -
96 Sqrt[2] E^(4/TeV) π^6 TeV^6) #1^3 + (-32 E^(
2/TeV) π^3 TeV^3 + 385 E^(4/TeV) π^6 TeV^6 +
32 Sqrt[2] E^(4/TeV) π^6 TeV^6) #1^4 + (64 E^(
2/TeV) π^3 TeV^3 -
512 Sqrt[2] E^(2/TeV) π^3 TeV^3) #1^5 + (-32 E^(
2/TeV) π^3 TeV^3 + 512 Sqrt[2] E^(2/TeV) π^3 TeV^3) #1^6 +
256 #1^8 &, 1],
y -> Root[-2 E^(
2/TeV) π^3 TeV^3 Root[
E^(4/TeV) π^6 TeV^6 -
64 E^(6/TeV) π^9 TeV^9 + (-4 E^(4/TeV) π^6 TeV^6 -
32 Sqrt[2] E^(4/TeV) π^6 TeV^6 +
128 E^(6/TeV) π^9 TeV^9) #1 + (390 E^(4/TeV) π^6 TeV^6 +
96 Sqrt[2] E^(4/TeV) π^6 TeV^6 -
64 E^(6/TeV) π^9 TeV^9) #1^2 + (-772 E^(
4/TeV) π^6 TeV^6 -
96 Sqrt[2] E^(4/TeV) π^6 TeV^6) #1^3 + (-32 E^(
2/TeV) π^3 TeV^3 + 385 E^(4/TeV) π^6 TeV^6 +
32 Sqrt[2] E^(4/TeV) π^6 TeV^6) #1^4 + (64 E^(
2/TeV) π^3 TeV^3 -
512 Sqrt[2] E^(2/TeV) π^3 TeV^3) #1^5 + (-32 E^(
2/TeV) π^3 TeV^3 +
512 Sqrt[2] E^(2/TeV) π^3 TeV^3) #1^6 + 256 #1^8 &, 1]^2 +
4 E^(2/TeV) π^3 TeV^3 Root[
E^(4/TeV) π^6 TeV^6 -
64 E^(6/TeV) π^9 TeV^9 + (-4 E^(4/TeV) π^6 TeV^6 -
32 Sqrt[2] E^(4/TeV) π^6 TeV^6 +
128 E^(6/TeV) π^9 TeV^9) #1 + (390 E^(4/TeV) π^6 TeV^6 +
96 Sqrt[2] E^(4/TeV) π^6 TeV^6 -
64 E^(6/TeV) π^9 TeV^9) #1^2 + (-772 E^(
4/TeV) π^6 TeV^6 -
96 Sqrt[2] E^(4/TeV) π^6 TeV^6) #1^3 + (-32 E^(
2/TeV) π^3 TeV^3 + 385 E^(4/TeV) π^6 TeV^6 +
32 Sqrt[2] E^(4/TeV) π^6 TeV^6) #1^4 + (64 E^(
2/TeV) π^3 TeV^3 -
512 Sqrt[2] E^(2/TeV) π^3 TeV^3) #1^5 + (-32 E^(
2/TeV) π^3 TeV^3 +
512 Sqrt[2] E^(2/TeV) π^3 TeV^3) #1^6 + 256 #1^8 &, 1] #1 -
2 E^(2/TeV) π^3 TeV^3 #1^2 + #1^4 &, 2]}, {x ->
Root[E^(4/TeV) π^6 TeV^6 -
64 E^(6/TeV) π^9 TeV^9 + (-4 E^(4/TeV) π^6 TeV^6 -
32 Sqrt[2] E^(4/TeV) π^6 TeV^6 +
128 E^(6/TeV) π^9 TeV^9) #1 + (390 E^(4/TeV) π^6 TeV^6 +
96 Sqrt[2] E^(4/TeV) π^6 TeV^6 -
64 E^(6/TeV) π^9 TeV^9) #1^2 + (-772 E^(4/TeV) π^6 TeV^6 -
96 Sqrt[2] E^(4/TeV) π^6 TeV^6) #1^3 + (-32 E^(
2/TeV) π^3 TeV^3 + 385 E^(4/TeV) π^6 TeV^6 +
32 Sqrt[2] E^(4/TeV) π^6 TeV^6) #1^4 + (64 E^(
2/TeV) π^3 TeV^3 -
512 Sqrt[2] E^(2/TeV) π^3 TeV^3) #1^5 + (-32 E^(
2/TeV) π^3 TeV^3 + 512 Sqrt[2] E^(2/TeV) π^3 TeV^3) #1^6 +
256 #1^8 &, 2],
y -> Root[-2 E^(
2/TeV) π^3 TeV^3 Root[
E^(4/TeV) π^6 TeV^6 -
64 E^(6/TeV) π^9 TeV^9 + (-4 E^(4/TeV) π^6 TeV^6 -
32 Sqrt[2] E^(4/TeV) π^6 TeV^6 +
128 E^(6/TeV) π^9 TeV^9) #1 + (390 E^(4/TeV) π^6 TeV^6 +
96 Sqrt[2] E^(4/TeV) π^6 TeV^6 -
64 E^(6/TeV) π^9 TeV^9) #1^2 + (-772 E^(
4/TeV) π^6 TeV^6 -
96 Sqrt[2] E^(4/TeV) π^6 TeV^6) #1^3 + (-32 E^(
2/TeV) π^3 TeV^3 + 385 E^(4/TeV) π^6 TeV^6 +
32 Sqrt[2] E^(4/TeV) π^6 TeV^6) #1^4 + (64 E^(
2/TeV) π^3 TeV^3 -
512 Sqrt[2] E^(2/TeV) π^3 TeV^3) #1^5 + (-32 E^(
2/TeV) π^3 TeV^3 +
512 Sqrt[2] E^(2/TeV) π^3 TeV^3) #1^6 + 256 #1^8 &, 2]^2 +
4 E^(2/TeV) π^3 TeV^3 Root[
E^(4/TeV) π^6 TeV^6 -
64 E^(6/TeV) π^9 TeV^9 + (-4 E^(4/TeV) π^6 TeV^6 -
32 Sqrt[2] E^(4/TeV) π^6 TeV^6 +
128 E^(6/TeV) π^9 TeV^9) #1 + (390 E^(4/TeV) π^6 TeV^6 +
96 Sqrt[2] E^(4/TeV) π^6 TeV^6 -
64 E^(6/TeV) π^9 TeV^9) #1^2 + (-772 E^(
4/TeV) π^6 TeV^6 -
96 Sqrt[2] E^(4/TeV) π^6 TeV^6) #1^3 + (-32 E^(
2/TeV) π^3 TeV^3 + 385 E^(4/TeV) π^6 TeV^6 +

32 Sqrt[2] E^(4/TeV) π^6 TeV^6) #1^4 + (64 E^(
2/TeV) π^3 TeV^3 -
512 Sqrt[2] E^(2/TeV) π^3 TeV^3) #1^5 + (-32 E^(
2/TeV) π^3 TeV^3 +
512 Sqrt[2] E^(2/TeV) π^3 TeV^3) #1^6 + 256 #1^8 &, 2] #1 -
2 E^(2/TeV) π^3 TeV^3 #1^2 + #1^4 &, 2]}} *)


The Plots for the two solution sets are

Column[
Plot[Evaluate[{x, y} /. sol[[#[[1]]]]], {TeV, 0.1, 10},
PlotRange -> #[[2]],
Frame -> True,
FrameLabel -> (Style[#, 14, Bold] & /@ {"TeV", ""}),
ImageSize -> 360,
PlotLegends -> Placed[{x, y}, {.5, .25}]] & /@
{{1, Automatic}, {2, All}}]
`

Your actual parameter values and the context of your problem would presumably tell you which solution set(s) is(are) of interest.

• Thank you very much for your help! I confirm it works well now. So my problem was mostly linked to the correct syntax to use (as I suspected). If I move from NSolve to Solve, how should I modify your code? Thanks again! Commented Oct 10, 2019 at 18:22