How to optimize code for obtaining large data set?

In an application, I need to find a region where some conditions must be true. However, these conditions are made of interpolating function solved by ParametricNDSolve. I will give a toy model example to clarify

eq1 = y1'[x]==y1^2+y2^4
eq2 = y2'[x]==2*y1^4-y2^2
soln=ParametricNDSolve[{eq1,eq2,y1[0]==a,y2[0]==b},{y1,y2},{x,0,1},{a,b}];
y1sln[a_,b_,x_]:=Evaluate[y1[a,b][x]/.soln];
y2sln[a_,b_,x_]:=Evaluate[y2[a,b][x]/.soln];


These run super fast and I encounter no problems at all. After that, there are the conditions, which I write as

c1[a_,b_,x_]:=y1sln[a,b,x]^2+y2sln[a,b,x];
c2[a_,b_,x_]:=y2sln[a,b,x]+y1sln[a,b,x];


Supposing that I need thoses conditions to be less or equal than zero, I want to find the region where the {a,b} values are allowed and plot them. To do that, I'm using the following

Conditions[a_, b_, x_] := Apply[And,Table[ToExpression["c" <> ToString[i] <> "[a,b,x]" <> "\[LessEqual]0"], {i, 1,
2}]]
list = Reap[Table[Sow[Conditions[a, b, x]], {a, 0.,1., 0.02}, {b, 0., 1., 0.02]]


Essentially this is the code I'm using to generate the list of points that I need to plot. However, it's taking a huge amount of time. I would like to know if there is a proper and faster way to do this.

• Your example doesn't work for me in its current form. As for the performance aspect: Keep in mind that ParametricNDSolve doesn't really solve anything when you evaluate the first block of code. It just generates the ParametricFunction expressions that do the actual solving, which is why the last step is so slow. Commented Aug 9, 2022 at 12:47
• Any advice to overcome using ParametricNDSolve or ParametricFunction? Commented Aug 9, 2022 at 12:54

OP seems to be open to different approaches, so I thought I try a graphical approach. Here I plot the vector field of the autonomous differential equation, and the region where both c1 and c2 are not positive:

{y1min,y1max}={-2,2};
{y2min,y2max}={-2,2};
Show[{
VectorPlot[{y1^2+y2^4,2*y1^4-y2^2},{y1,y1min,y1max},{y2,y2min,y2max},
VectorScaling->Automatic,VectorPoints->Fine,
FrameLabel->{Subscript["y",1],Subscript["y",2]}],
RegionPlot[And[y1^2+y2<=0,y2+y1<=0],{y1,y1min,y1max},{y2,y2min,y2max}],
Graphics[{Red,Dashing[0.03],Line[{{y1min,0},{y1max,0}}]}]
}]


I also plot a dashed red line that I will come to in a moment. Here is the plot:

OP wants to know for which initial points $$y_1(0) = a$$, $$y_2(0) = b$$ the solution $$y_1(x)$$, $$y_2(x)$$ is in the blue region. Unfortunately, OP did not specify for which $$x$$ this conditions must hold. For all $$x$$? For at least one $$x$$?

We can nevertheless say something conclusive if we restrict (as OP did in the code in the question) to the case $$b \geq 0$$, which means that we start somewhere above the red dashed line, and to $$x \geq 0$$. In that case we never go below the red dashed line for $$x \geq 0$$, in particular we never enter the blue region. This can be seen from the plot, since along the red dashed line, the arrows always point upward. This can also be seen from the differential equation itself, by setting $$y_2=0$$ there. (The only exception is $$a=b=0$$, where the solution remains at the origin eternally, and which belongs to the blue region.)

Comment. Note that solutions do not necessarily exist for all $$x$$. It can happen that the solution escapes to infinity after finite time $$x$$, because the differential equation is nonlinear.

• Thank you, it was clarifying Commented Aug 11, 2022 at 0:52