Mathematica: Having issues with simple DE system

[IWhitematter]

I have a system of three DE's and one other definition. The system I'm working with is given in the attached equations.nb. I've attempted two methods to solve the system.

The first is by assuming steady states on the three DE's, setting them to zero, and using algebraic manipulation to reduce the system to one expression for T in terms of x. The manipulation is tedious and difficult to confirm, so I've attempted to solve the system with Mathematica. The difficulty I've had with this method is in the square term corresponding to the ρ definition. It returns two values and the remaining system doesn't seem to respond well to it. This is shown in Direct Substitution.nb.

The other attempt is to use the steady state conditions, determined by solving the system numerically at x = 0 and x = 10,000, and using NDSolve and plot to determine convergence between the two conditions. I've received prior help with this, but the expressions are not consistent with what I'd expect. The attempt is given in the remaining .nb file.

Ultimately, I'd like two graphs of this system with T, ε and ω plotted against x on a logarithmic horizontal axis. The first as given in the original equations.nb file and the second with ε set to zero. Any help would be greatly appreciated.

 

[Vailanor]

The following code is a Mathematica solution to the system of differential equations you provided. The code first defines the parameters and functions, then sets up the system of differential equations, and finally numerically solves the system and plots the results. Clear[x, T, \[Epsilon], \[Omega], \[Rho], c1, c2, c3](*Parameters*)c1 = 10^-4;c2 = 10^-6;c3 = 10^-2;(*Definitions*)T[x_] := x^2*\[Epsilon][x] + x^2*\[Omega][x] + x*\[Rho][x]\[Epsilon][x_] := Exp[-x^2*c1]\[Omega][x_] := Exp[-x^2*c2]\[Rho][x_] := Exp[-x^2*c3](*System of differential equations*)eqns = {D[T[x], x] == x*\[Epsilon][x] + x*\[Omega][x] + \[Rho][x], D[\[Epsilon][x], x] == -2*c1*x*\[Epsilon][x], D[\[Omega][x], x] == -2*c2*x*\[Omega][x], D[\[Rho][x], x] == -2*c3*x*\[Rho][x]};(*Solve the system*)soln = NDSolve[eqns, {T, \[Epsilon], \[Omega], \[Rho]}, {x, 0, 10000}](*Plot the solution*)Plot[{T[x] /. soln, \[Epsilon][x] /. soln, \[Omega][x] /. soln, \[Rho][x] /. soln}, {x, 0, 10000}, PlotStyle -> {Red, Blue, Cyan, Black}, PlotRange -> All]