# Solving a system of coupled Nonlinear ODEs using numeric and get ndnum error

I would like to reproduce the figures 2 of section 4.1 of https://arxiv.org/pdf/1401.6088.pdf. I trying to solve a system of coupled ODEs (4.4 and 4.5) and then plotted them . The system is $$\frac{\partial \left(v'(x) z(x)^{2 (1-\xi )} F(v,z)+z(x)^{1-\xi } z'(x)\right)}{\partial x}-\frac{1}{2} v'(x)^2 z(x)^{2 (1-\xi )} \frac{\partial F(v,z)}{\partial v}=0$$ $$-\frac{\text{d\theta } \left(-v'(x)^2 z(x)^{2 (1-\xi )} F(v,z)-2 v'(x) z(x)^{1-\xi } z'(x)+1\right)}{z(x)}-(1-\xi ) v'(x) z(x)^{-\xi } \left(v'(x) z(x)^{1-\xi } F(v,z)+z'(x)\right)-\frac{1}{2} v'(x)^2 z(x)^{2 (1-\xi )} \frac{\partial F(v,z)}{\partial z}+\frac{\partial \left(v'(x) z(x)^{1-\xi }\right)}{\partial x}=0$$ where the prime denotes derivative w.r.t x and the constants and functions are defined as $$\text{d\theta }=d-\theta ;\theta =\frac{2}{3};\xi =1.5;d=2;a=0.5;l=8;\epsilon =\frac{1}{10^{16}};$$ $$F(\text{v\_},\text{z\_})\text{:=}1-M(v) z(x)^{\text{d\theta }+\xi }$$ $$M(\text{v\_})\text{:=}\frac{1}{2} \left(\tanh \left(\frac{v(x)}{a}\right)+1\right)$$ The boundary conditions are $$v\left(\frac{l}{2}\right)=-20,z\left(\frac{l}{2}\right)=\frac{1}{10^3},z(\epsilon )=5,z'(\epsilon )=0,v'(\epsilon)=0$$ The Codes are

 eq1=D[z[x]^(1 - \[Xi]) (z[x]^(1 - \[Xi]) F v'[x] + z'[x]), x] - z[x]^(2 (1 - \[Xi])) D[F, v[x]] v'[x]^2/2 == 0;

eq2= D[z[x]^(1 - \[Xi]) v'[x], x] -d\[Theta]/
z[x] (1 - F z[x]^(2 (1 - \[Xi])) v'[x]^2 -
2 z[x]^(1 - \[Xi]) z'[x] v'[x]) -
z[x]^(2 (1 - \[Xi])) D[F, z[x]] v'[x]^2/
2 - (1 - \[Xi]) z[x]^-\[Xi] (z'[x] + z[x]^(1 - \[Xi]) F v'[x]) v'[
x] == 0;

d\[Theta]=d-\[Theta];\[Theta]=2/3;\[Xi]=1.5;d=2;a=0.5;l=8;\[Epsilon]=\$MachineEpsilon;
bcs={v[l/2]==-20,z[l/2]==10^-3,z[\[Epsilon]]==5,z'[\[Epsilon]]==10^-9};

F[v_,z_]:=1-M[v]z[x]^(d\[Theta]+\[Xi]);
M[v_]:=1/2 (1+Tanh[v[x]/a]);


I trying to solve by following command but we get many errors including infinite expressions and ndnum error

NDSolve[{eq1,eq2,bcs},{z,v},{x,\[Epsilon],4}]


NDSolve::ndnum: Encountered non-numerical value for a derivative at x == 2.220446049250313*^-16.
I don't have any idea how to solve this system. Please give me some hint.

Thanks

• have you seen NDSolve ? Perhaps obvious..but since you don't say what you tried its hard to see where you are stuck. – george2079 May 16 '17 at 15:18
• what is M on the right side of the expression defining the function M` ? – george2079 May 16 '17 at 15:22
• Check out mathematica.stackexchange.com/questions/146049/… for setting StartingInitialConditions when trying to solve a boundary value problem. I can't immediately find a set of initial values which does let Mathematica fit your boundary conditions though, perhaps you have more idea what they should start at. – KraZug May 17 '17 at 12:16
• Other issues: You haven't specified t in this code here, and your definition for M probably wants to be a pattern (using an underscore) so it works for general v, not just exactly v[x]. And your range goes to 5, not just l/2. And you introduce decimals via [Xi] and a, it is probably better to use fractions here. – KraZug May 17 '17 at 12:17