# Problem with using NDSolve to solve coupled differential equations [closed]

I am having a problem with the following code:

rmax = 100, g = 3, n = 10^10, B = 160^4, T = 170, r0 = 1.477, L = n/(4/3*pi*r^3)
a = 1/(20*r0^3), mu = ((6*L*pi^2)/g)^1/3

Er(r) = 3*pr(r) + 4*B - 16/27*mu*T(mu^2 - 2*mu*T + T^2)

s =
NDSolve[
{m'[r] == a*r^2 Er[r],
pr'[r] == -(r0/r)*((pr[r] + Er[r]) (m[r] + a*r^3* pr[r]))/(r - 2 r0*m[r]),
m[1] == 1, pr[1] == 1},
{m, pr}, {r, 1, rmax}]


It is showing the error

NDSolve : Encountered non-numerical value for a derivative at r == 1.

How to solve the problem raised by the message?

## closed as off-topic by bbgodfrey, user21, MarcoB, Alex Trounev, mikuszefskiJun 21 at 8:13

This question appears to be off-topic. The users who voted to close gave this specific reason:

• "This question arises due to a simple mistake such as a trivial syntax error, incorrect capitalization, spelling mistake, or other typographical error and is unlikely to help any future visitors, or else it is easily found in the documentation." – bbgodfrey, user21, MarcoB, Alex Trounev, mikuszefski
If this question can be reworded to fit the rules in the help center, please edit the question.

• There are many errors in your code. To begin with, use ";" instead of "," to separate assignements. Use Pi instead of pi. Make sure functions are defined with "[ ]" instead of "( )". – Whelp Jun 18 at 8:09
• Check your code: Er(r)->Er[r], pr(r)->pr[r], pi->Pi. The corrected system equations show an infinite slope at r==1, that's why NDSolve doesn't find a solution! – Ulrich Neumann Jun 18 at 8:10
• Is it a stiff equation kind of thing? Can I plot the function any how??? – Sovan Jun 19 at 6:49

### Edit

This is not an answer, but an extended comment which reveals the syntax errors that must be corrected before progress can be made on the real problem inherent in the OP's system of equations.

With correct Mathematica syntax, there is no complaint about non-numerical quantities but there is still an error.

rmax = 100; g = 3; n = 10^10; B =  160^4; T = 170 ; r0 = 1.477; L = n/(4/3 Pi r^3);
a = 1/(20 r0^3); mu = ((6 L Pi^2)/g)^1/3;

Er[r_] = 3 pr[r] + 4 B - 16/27 mu T (mu^2 - 2 mu T + T^2);

{mF, prF} =
NDSolveValue[
{m'[r] == a r^2 Er[r],
pr'[r] == -(r0/r) ((pr[r] + Er[r]) (m[r] + a r^3 pr[r]))/(r - 2 r0 m[r]),
m[1] == 1, pr[1] == 1},
{m, pr}, {r, 1, rmax}];


NDSolveValue::ndsz: At r == 1., step size is effectively zero; singularity or stiff system suspected.

### Notes

1. pi should be Pi.
2. Multiple assignments should be separated with ; not ,.
3. E(r) needs to be replaced with E[r_].
4. I recommend NDSolveValue over NDSolve because it returns a list of interpolation functions. I find that result easier to use in further calculations.
• You have missed the L in mu expression while writing the code. – Sovan Jun 18 at 21:03
• After inserting mu expression with the missing L it is showing the same problem. – Sovan Jun 19 at 6:50
• @Sovan. I missed both the L in the definition of mu and the pi in the definition of L. I have made the required edits and revised my conclusions accordingly. – m_goldberg Jun 19 at 15:06
• To address the "At r==1, step size is effectively zero..." issue, the first thing I'd try is adding WhenEvent[r==1,"CrossSlidingDiscontinuity"] and see if that fixes it. – Matt Jun 19 at 17:25

Using code @m_goldberg and explicit Euler, we can get a solution.

rmax = 100; g = 3; n = 10^10; B = 160^4; T = 170; r0 = 1477/100; L =
n/(4/3 Pi r^3);
a = 1/(20 r0^3); mu = ((6 L Pi^2)/g)^1/3;

Er[r_] = 3 pr[r] + 4 B - 16/27 mu T (mu^2 - 2 mu T + T^2);

{mF, prF} =
NDSolveValue[{m'[r] == a r^2 Er[r],
pr'[r] == -(r0/r) ((pr[r] + Er[r]) (m[r] + a r^3 pr[r]))/(r -
2 r0 m[r]), m[1] == 1, pr[1] == 1}, {m, pr}, {r, 1, rmax},
StartingStepSize -> 10^-2, WorkingPrecision -> 30,
Method -> "ExplicitEuler"];

{Plot[mF[r], {r, 1, rmax}], Plot[prF[r], {r, 1, rmax}]}
`