# How can i solve and plot this non linear ODEs usin NDSOLVE?

I am trying to solve the following non-linear ODEs

Eqn1 = n*(-f''[x])^(n - 1)f'''[x] - mf'[x]^2 + (m*(2*n-1)+1)/(n+1)f[x]f''[x]+ M2^2f'[x] == 0; ,
Eqn2 = n(-theta'[x])^(n - 1)theta''[x] + (m(2 *n - 1) + 1)/(n + 1)f[x]theta'[x] - mf'[x]theta[x] == 0;


where m,n, alpha, b and M2 are parametes, with 5 boundary conditions

f[0] == 0, f'[0] == 1 + alpha f''[0], f'[N1] == 0, theta[0] == 1 + b theta'[0], theta[N1] == 0


• Please enter the code in a more readable format. You can read about the formatting options here: Markdown help#Code and Preformatted Text Commented Sep 11, 2020 at 9:15
• What does "  BC3 = f'[N1] == 0 " mean? Further, assignments are written with 1 and equation with 2 equal signs. Commented Sep 11, 2020 at 9:41

Ok, so there are a lot of quantities in the OP that are not defined and/or discussed but I used some made up values to demonstrate the procedure.

I am changing theta to g for my own convenience.

Set up the equations

Eqn1 = n*(-f''[x])^(n - 1) f'''[x] -
m f'[x]^2 + (m*(2*n - 1) + 1)/(n + 1) f[x] f''[x] + M2^2 f'[x] ==
0; Eqn2 =
n (-g'[x])^(n - 1) g''[x] + (m (2*n - 1) + 1)/(n + 1) f[x] g'[x] -
m f'[x] g[x] == 0;


Setting up the made up values for the parameters

m = 1;
n = 2;
alpha = 1;
b = 1;
M2 = 2;
N1 = 0;


You solve Eqn1 with the following piece of code

sltn = NDSolve[{Eqn1, f[0] == 0, f'[0] == 1 + alpha f''[0],
f'[N1] == 0}, f[x], {x, 0, 10},
Method -> {"StiffnessSwitching", "NonstiffTest" -> False,
Method -> {"ExplicitRungeKutta", Automatic}}, AccuracyGoal -> 5,
PrecisionGoal -> 5, MaxSteps -> Infinity] // Flatten


And you can actually plot the solution to look at it

Plot[f[x] /. sltn, {x, 0, 10}]


Of course, without further clarification, not much progress can be made but hopefully the above is a basic guiding principle.

Edit: Let me demonstrate the solution for both functions. I am doing this for clarity.

With the above initializations you run the following

sltn = NDSolve[{Eqn1, Eqn2, f[0] == 0, f'[0] == 1 + alpha f''[0],
f'[N1] == 0, g[N1] == 0, g[0] == 1 + b g'[0]}, {f[x], g[x]}, {x,
0, 10}, Method -> {"StiffnessSwitching", "NonstiffTest" -> False,
Method -> {"ExplicitRungeKutta", Automatic}}, AccuracyGoal -> 5,
PrecisionGoal -> 5, MaxSteps -> Infinity] // Flatten


And then you can plot your solutions

Plot[f[x] /. sltn[[1]], {x, 0, 10}, PlotRange -> {{0, 10}, {-1, 1}}]
Plot[g[x] /. sltn[[2]], {x, 0, 10}, PlotRange -> {{0, 10}, {-10, 1}}]


If you want to solve for more general values of the parameters you might want to look up the ParametricNDSolve command here

Hope this helps a bit.

• Probably N1>0 is upper limit of x. sltn = NDSolve[......, {x, 0, N1}.... Commented Sep 11, 2020 at 11:45
• I do not disagree with you, however, there are not enough clarifications in the OP and my intention was to provide some code with the basic guiding principles so the user can apply it to the actual problem
– user49048
Commented Sep 11, 2020 at 11:47