NDSolve Plotting issue

Question

I am trying to solve a system of ODEs with one extra boundary condition.

eq1 = f'''[y] + f[y] f''[y] - m f'[y] - f'[y]^2 - s (f'[y] + y/2 f''[y]) == 0;
eq2 = g''[y] - (s/2(3 g[y] + y g'[y]) + 2 g[y] f'[y] - g'[y] f[y]) + f''[y]^2 == 0;

With boundary conditions

bc = {f[0] == 0, f'[0] == 1, g[0] == 1, f''[b] == 0, g'[b] == 0};

With the arbitrary trial values, the integrated solutions will generally not satisfy the outer boundary conditions.

f''[b] == 0 and f[b] == s b/2

I need to find a way to adjust the two trial values such that the numerical solution eventually matches the required boundary conditions and satisfy the relation f[b] = s b/2.

Here are my questions

• If I fix m (say = 1), then how I can plot b vs s?
• If I fix s (say = 1.2 and 0.8), then how I can plot f'[b] vs m?
• If I fix s (say = 1.2 and 0.8), then how I can plot f''[0] vs m?

What value b will take?

Here is my try, but it's not working

g[b_?NumericQ, m_?NumericQ, s_?NumericQ] := 
  f[b] /. 
    NDSolve[{eq1, eq2, f[0] == 0, f'[0] == 1, g[0] == 1, f''[b] == 0, g'[b] == 0}, 
      {f, g}, {y, 0, b}]

s2[b_, s_] := FindRoot[g[b, s, m], {b, 0}]

Paper and Pen implementation of Shooting method As it is mentioned by @Peltio, here is the possible Shooting method.

Note $theta$ is g and $\eta$ is y.

The output should look like this

1. b vs s

2. f'[b] vs m

3. f''[0] vs m

Answer 1

The following does work in principle but gives a lot of warnings and I don't think you should trust the numbers it gives, although they roughly seem to reproduce the first plot (which is the only I checked but don't know the m it was produced with). I have tried to keep this simple but robust and try to show some techniques that I found helpful in the past. So I hope it can serve as a guide about how to access similar problems.

Avoid Redefinitions

to keep things simple I just do a

ClearAll[f, θ, η, f0, f1, f2, θ0, θ1, s, m]

in the beginning. If you plan to pack all that into a function I would certainly make all those local variables of a scoping construct, e.g. Module...

Derive Equations

One of the probably underrated features of Mathematica is that you can input in almost exactly the form you find it in textbooks or publications. It will help a lot to avoid unnecessary typos if one makes good use of that. Thus I even sticked to the greek letters (which now even look nice here and can be copied to Mathematica thanks to shrx editing). To check with the original work you want to convert the results to StandardForm within Mathematica, e.g. via Cell->Convert To->StandardForm. It probably is worth noting that I detected such a type when trying this out before posting this, that typo is now fixed (missing (f1[η])^2 in first equation). I also let Mathematica do the simple algebra, which helps to avoid unneccessary typing and many of those annoying errors that people mentioned in the comments to your question. Here I was using this to define the equations:

eqns[m_, s_] = Function[
  {f, θ, η},
  Evaluate[{
     f2'[η] == s*(f1[η] + η/2 f2[η]) ] + (f1[η])^2 - f0[η]*f2[η] + m*f1[η],
     θ1'[η] == (s/2 (3*θ0[η] + η*θ1[η]) + 2*θ0[η] f1[η] - θ1[η]*f0[η]) - (f2[η])^2
     } /. {
     f0 -> f, f1 -> Derivative[1][f], f2 -> Derivative[2][f],
     θ0 -> θ, θ1 -> Derivative[1][θ]
     }
   ]
  ]

bc[β_] = Function[{f, θ, η},
  Evaluate[{
    f[0] == 0, f'[0] == 1, θ[0] == 1,
    f''[β] == 0, θ'[β] == 0
    }]
  ]

note that by using Derivative one can write those replacments in a very elegant and robust way. I also consider it good practive to make the parameters function arguments in such cases.

Find Solution

I split this in two parts, first is to solve the boundary condition problem. In newer versions of Mathematica NDSolve has a shooting method built in and will choose that automatically if given according boundary conditions. Thus this will give you a solution of f in the form of an InterpolatingFunction:

fSolution[β_?NumericQ, m_?NumericQ, s_?NumericQ] := 
 First[NDSolveValue[
   Evaluate[Flatten[{eqns[m, s][f, θ, η], bc[β][f, θ, η]}]],
   {f}, {η, 0, β}
 ]]

This is a good stage where to test whether everything works as desired, e.g. by evaluating:

fSolution[1, 1, 1]

you would typically also investigate whether that solution fullfills the boundary conditions and in general looks like an appropriate solution...

Having this, we can now try to find that β which will fullfill the additional condition:

βSolution[m_?NumericQ, s_?NumericQ] := Module[{β},
  β /. First[FindRoot[fSolution[β, m, s][β] == (s*β)/2, {β, 0.01, 10}]]
  ]

the initial values of FindRoot here are adjusted for the following table/plot, you probably might need to do something smarter to make this work for aribtrary cases...

Calculate and Plot values

As the above is using two nested FindRoot (one with two variables implicitly by NDSolves shooting method) and an NDSolve within that calculating a solution isn't very fast. In such cases it often makes much more sense to decide at which points one wants to evaluate and calculate them instead of letting Plot do it's magic which in general will lead to a lot more evaluations. So here is quite rough approximation which is mainly good to check whether the above does what it should, for a nice plot you'd want to use some more points:

βOverS = Module[{x},
  Table[x = {s, βSolution[1, s]}; Print[x]; x, {s, 0.2, 1, 0.1}]
  ]

ListLinePlot[βOverS, PlotRange -> {0, 10}, Frame -> True]

I have introduced a very crude monitoring of the progress. As mentioned above both FindRoot and NDSolve give some warnings about failing to find a good result and it will not even work for smaller values of s. I guess that you could get better results and get them faster with adjusting the various methods of the two functions. But that is a wide field and probably would justify a new question which should only concentrate on that aspect...

EDIT: additional plots

first, to make the following more comfortable one might want to memoize all beta-Values that were calculated. This can be done by redefining βSolution like this:

βSolution[m_?NumericQ, s_?NumericQ] := βSolution[m,s] = Module[{β},
  Print[βSolution->{m,s}];
  β /. First[FindRoot[fSolution[β, m, s][β] == (s*β)/2, {β, 0.01, 10}]]
]

I also added a Print so one can immediately see what is currently calculated. With these preparations one can then (almost) reproduce the other two plots like this:

fPrimeOverM[s_] := Module[{β},
    Table[
        β = βSolution[m, s];
        {m, Derivative[1][fSolution[β, m, s]][β]},
        {m, 0, 8, 0.5}
    ]
];

ListLinePlot[{fPrimeOverM[1.2], fPrimeOverM[0.8]},
    Frame -> True,PlotStyle -> {Red, Black},
    PlotRange -> {0, 0.5},AspectRatio -> 0.8
]   

and:

fSecondOverM[s_] := Module[{β},
  Table[
   β = βSolution[m, s];
   {m, -Derivative[2][fSolution[β, m, s]][0]},
   {m, 0, 8, 0.5}
   ]
  ]

ListLinePlot[Evaluate[{fSecondOverM[1.2], fSecondOverM[0.8]}],
 Frame -> True,
 PlotStyle -> {Red, Black},
 PlotRange -> {0, 3.5},
 AspectRatio -> 0.8]

On my machine that gives larger discrepancies for small values of M, and that is also where FindRoot and NDSolve give a lot of warnings which indicates that these should be taken seriously. To really get reliable results you would definitely have to improve the settings of especially NDSolve (and of course recheck all definitions, including potential typos in the original work).

Answer 2

This is not the answer to your problem, but perhaps this can help you find a solution to your plotting problem. I would start with defining the equations in terms of the parameters s and m:

eq1[m_, s_] = f'''[y] + f[y] f''[y] - m f'[y] - f'[y]^2 - s (f'[y] + y/2 f''[y]) == 0;
eq2[m_, s_] = g''[y] - (s/2(3 g[y] + y g'[y]) + 2 g[y] f'[y] - g'[y] f[y]) + f''[y]^2 == 0;

Then I will solve the systems of ODEs with the condition at one extreme of the interval.

vals[m_, s_, b_] := (
    {f1, g1} = ({f, g} /. NDSolve[{eq1[m, s], eq2[m, s],
                f[0] == 0, f'[0] == 1, f''[0] == -2,
                g[0] == 1, g'[0] == -.2
                }, {f, g}, {y, 0, b}][[1]]);
    {f1[b], g1[b]})
(*You can choose to return the derivatives if you wish*)

Mind you: this is not your system of ODEs+BCs. I had to make up further initial conditions for f and g. Your problem is more complicated in that you need to impose conditions at three points. (If m and s were known, this would be a straightforward problem of shooting to a fitting point - you would solve it using the extremes as initial/final condition and use the solution to compute the value at the same common point).

Now, if you specify one of the three values, you can plot the solutions in terms of the remaining twos. For example, suppose the middle point is b = 2, then the values of f and g in 2 can be seen as functions of m and s:

{Plot3D[vals[m, s, 2][[1]], {m, -3, 3}, {s, -3, 3}], 
    Plot3D[vals[m, s, 2][[2]], {m, -3, 3}, {s, -3, 3}]}

Incidentally, if you can get the solution to your system, this will allow you to do what you ask. Just specify the fixed value and plot vals[[1]] - or just vals if you choose to return a single function's value - versus the other two. Of course vals can be made to return f, f' or f''.

I hope this can help you getting started.