How to remove the warnings and errors from this manipulate code?

Question

I'm experiencing several problems with the manipulate code below. This code numerically solves three first-order non-linear differential equations, and then output a 3D plot of the curve, using the initial conditions entered by the user from the three sliders :

Clear["Global`*"]
Theta1 = -20;
Theta2 = 20;

dynamics[om_, or_, ol_] := dynamics[om, or, ol] = NDSolve[{
    oM'[theta] == (oM[theta] + 2oR[theta] - 2oL[theta] - 1)oM[theta],
    oR'[theta] == (oM[theta] + 2oR[theta] - 2oL[theta] - 2)oR[theta],
    oL'[theta] == (oM[theta] + 2oR[theta] - 2oL[theta] + 2)oL[theta],

    oM[0] == om, oR[0] == or, oL[0] == ol

    }, {oM, oR, oL}, {theta, Theta1, Theta2},
    Method -> "StiffnessSwitching"
    (* constraints to be imposed : oM, oR, oL > 0 only *)
]

Tmin[om_, or_, ol_] := Theta1 (* to be fixed *)
Tmax[om_, or_, ol_] := Theta2 (* to be fixed *)

solution[om_, or_, ol_] := ParametricPlot3D[
    Evaluate[{oM[theta], oR[theta], oL[theta]}/.dynamics[om, or, ol]],
    {theta, Tmin[om, or, ol], Tmax[om, or, ol]}
]

Manipulate[
    Show[
    solution[om, or, ol],
    PlotRange -> {{0, 2}, {0, 2}, {0, 2}},
    SphericalRegion -> True
    ],
    {{om, 0.3, a}, 0, 2, 0.01},
    {{or, 0.0, b}, 0, 2, 0.01},
    {{ol, 0.7, c}, 0, 2, 0.01}
]

At compilation, I'm getting several error messages that I don't know how to solve. The manipulate box should be regular for all values entered from the sliders (from 0 up to 2 or more), especially when any of the variables is set to 0. The oM, oR and oL variables should only be positive, so a constraint should be added to the NDSolve part, and the curve extremities should be properly defined with the Tmin and Tmax definitions above. Currently, it doesn't work well.

How can I fix these problems ?

EDIT : Here are two typical messages that annoys me :

... step size is effectively zero; singularity or stiff system suspected

... lies outside the range of data in the interpolating function. Extrapolation will be used.

Answer 1

a bit of an extended comment, this is a minimal demonstration of the error you see. (You should have posted something like this to begin with) Note dynamic/manipulate has nothing to do with it, nor does the initial condition being in the middle of the domain cause the problem.

r = With[{om = 1, or = 2, ol = .3}, First@NDSolve[{
      oM'[ theta] == (oM[theta] + 2 oR[theta] - 2 oL[theta] - 1) oM[theta],
      oR'[theta] == (oM[theta] + 2 oR[theta] - 2 oL[theta] - 2) oR[theta],
      oL'[theta] == (oM[theta] + 2 oR[theta] - 2 oL[theta] + 2) oL[theta],
      oM[0] == om, oR[0] == or, oL[0] == ol},
     {oM, oR, oL},
     {theta, 0, 2}]];

NDSolve::ndsz: At theta == 0.35134263731177495`, step size is effectively zero; singularity or stiff system suspected. >>

Plot[{(oR /. r)[x], (oM /. r)[x], (oL /. r)[x]}, {x, 0, 1}, 
 PlotRange -> {0,40}]

the pde is simply blowing up.

increasing ol to .4 yields a stable solution:

Now that we see whats going on, go back to Manpulate, adding Quiet to surpress the error and extract the valid range from each solution:

Manipulate[
 Quiet[
  r = With[{om = 1, or = 2}, First@NDSolve[{
       oM'[theta] == (oM[theta] + 2 oR[theta] - 2 oL[theta] - 1) oM[
          theta],
       oR'[theta] == (oM[theta] + 2 oR[theta] - 2 oL[theta] - 2) oR[
          theta],
       oL'[theta] == (oM[theta] + 2 oR[theta] - 2 oL[theta] + 2) oL[
          theta],
       oM[0] == om, oR[0] == or, oL[0] == ol},
      {oM, oR, oL},
      {theta, 0, 2}]]];
  range = Flatten[(oM /. r)["Domain"]];
  Plot[{(oR /. r)[x], (oM /. r)[x], (oL /. r)[x]}, 
      Evaluate[Join[{x}, range]], PlotRange -> {0, 40}],
{{ol, .4}, .3, 1}]

Finally, fixing the original should look like this:

solution[om_, or_, ol_] := (
  result = Quiet@First@dynamics[om, or, ol];
  range = (oM /. result)["Domain"][[1]];
  ParametricPlot3D[{oM[theta], oR[theta], oL[theta]} /. result, 
   Evaluate[Join[{theta}, range]]])