# Numerically solving coupled ODE's with a parameter as initial condition

i'm currently trying to numerically solve a set of coupled ODE's to obtain the functions p(r), h(r) and m(r) in the range of r1 <= r <= r2 with initial conditions m(r1)=a=const and p(r1)=b=const. It is my goal to vary the parameters a and b to find the highest possible value for m(r2). So far i have no problems solving the ODE's for fixed a and b, but i have no idea how to vary the two parameters. I'm really new to Mathematica and after working around with ParametricNDSolve, ParametricPlot etc. for one week without making any progress i am now desperate for help.

Here's the code i made so far that solves the ODE's where i fixed r1=1.8, r2=2.2 and a=0.8, b=1. Please note that there is also another condition that a < r1/2, but i don't think that this is very important to my problem.

r1 = 1.8;
r2 = 2.2;
ode1 = {h'[r] == (1/r)*(1 - 8*Pi*p[r]*r^2 - h[r])};
ode2 = {p'[r] == -(2 p[r]/r) - (h'[r]/h[r])*p[r]*((1 + 8*Pi*p[r]*r^2 - 3*h[r])/(1 - 8*Pi*p[r]*r^2 - h[r]))};
{H, P, M} = {h, p, m} /.NDSolve[{ode1, ode2 , m[r] == (r/2)*(1 - h[r]), m[r1] == 0.8, p[r1] == 1}, {h, p, m}, {r, r1, r2}] // First


Thank you very much in advance.

• Welcome to Mathematica.SE! 1) As you receive help, try to give it too, by answering questions in your area of expertise. 2) Take the tour and check the faqs! 3) When you see good questions and answers, vote them up by clicking the gray triangles, because the credibility of the system is based on the reputation gained by users sharing their knowledge. Also, please remember to accept the answer, if any, that solves your problem, by clicking the checkmark sign!
– user9660
May 3, 2016 at 12:05

Oh! I'm coming to late!

r1 = 1.8;
r2 = 2.2;
ode1 = h'[r] == (1/r)*(1 - 8*Pi*p[r]*r^2 - h[r]);
ode2 = p'[
r] == -(2 p[r]/r) - (h'[r]/h[r])*
p[r]*((1 + 8*Pi*p[r]*r^2 - 3*h[r])/(1 - 8*Pi*p[r]*r^2 - h[r]));

sol = ParametricNDSolve[{ode1, ode2, m[r] == (r/2)*(1 - h[r]),
m[r1] == a, p[r1] == b}, {h, p, m}, {r, r1, r2}, {a, b}]


One way to find max m[a,b][r2]

MaximalBy[
Table[Evaluate[{a, b, m[a, b][2.2] /. sol}], {a, 0.2, 0.8, 0.1}, {b,
0, 0.5, 0.1}][[1]], Last]
(* {{0.2, 0.3, 0.970011}} *)

Manipulate[
Plot[Evaluate[{h[a, b][r], p[a, b][r], m[a, b][r]} /. sol], {r, r1,
r2}, PlotTheme -> "Detailed", PlotLegends -> {h, p, m}], {a, 0.1,
0.3, Appearance -> "Open"}, {b, 0.1, 0.5, Appearance -> "Open"}]


• Thank you very, very much for your answer! Don't worry, even though you answered later, your answer sure provides more useful information to me since you also gave me the means to find the maximum value! I really appreciate your help, you sure did me an enormous favor! Even though i need to do some research to understand what exactly you did, but reverse-engineering won't be too hard i guess :p
– Xeno
May 3, 2016 at 17:17
• @Xeno Of course, you can do it
– user36273
May 3, 2016 at 20:12
r1 = 1.8;
r2 = 2.2;
ode1 = {h'[r] == (1/r)*(1 - 8*Pi*p[r]*r^2 - h[r])};
ode2 = {p'[r] == -(2 p[r]/r) - (h'[r]/h[r])*
p[r]*((1 + 8*Pi*p[r]*r^2 - 3*h[r])/(1 - 8*Pi*p[r]*r^2 - h[r]))};

hpm = ParametricNDSolveValue[{ode1, ode2, m[r] == (r/2)*(1 - h[r]),
m[r1] == a, p[r1] == b}, {h, p, m}, {r, r1, r2}, {a, b}]

Plot[Evaluate[Through[hpm[0.8, 1.][t]]], {t, 1.8, 2.2}]


• Thank you very, very much for your help! Through, ParametricNDSolveValue? Those certainly are new commands to me! Thank you again for your help, i really appreciate it.
– Xeno
May 3, 2016 at 17:09