# How can I optimize plotting the solution to ParametricNDSolve as a function of a range of values of the parameter?

I'm using ParametricNDSolve to find a solution to a set of first order ODEs which have one independent variable, z. The set of initial values contains a function of t.

(*Constants*)
G = 0.00589667;
B = 2066.84;

ode = {y1'[z] == 2/3 * G * y2[z] * y3[z],
y2'[z] == -B * y3[z] - 2/3 * G * y1[z] * y3[z],
y3'[z] == B * y2[z]};

iv = {y1[0] == 0., y2[0] == 0., y3[0] == -Exp[-(t - 5)^2]};


ParametricNDSolve treats t as a parameter and returns three parametric interpolating functions with parameter t, independent variable z, and dependent variables y1, y2, and y3:

sol = ParametricNDSolve[{ode, iv}, {y1, y2, y3}, {z, 0., 50.}, t]

{y1 -> ParametricFunction[ <> ],
y2 -> ParametricFunction[ <> ],
y3 -> ParametricFunction[ <> ]}


Then I can evaluate these functions easily to obtain the y value at any t and z in range like this: y1[t][z] /. sol.

I want to plot the y values at a single value of z, over a range of t:

Plot[Evaluate@{y1[t][l] /. sol, y2[t][l] /. sol, y3[t][l] /. sol}, {t, 0, 10}, PlotRange -> All]


This does produce the plot I want, but it takes a very long time, over an hour.

I timed ParametricNDSolve[] and that is not the problem. I then tried a different method of plotting: I evaluated the functions for t and z first, then used ListPlot[]:

y1out = y1[#][50.] & /@ Range[0, 10, 0.01] /. sol;
ListPlot[y1out, PlotRange -> All]


But the improvement was negligible, so I believe it is the evaluation of the parametric functions that is inefficient.

I should note that using the Plot[] function is ideal because it changes its sampling rate with the slope of the curve. I can't predict the shape of the curves, and want to produce the shape of the curves reliably.

I am not a programmer and unfortunately don't have a Mathematica expert around to consult. I do have access to a high performance computer and so I can look into learning how to code in parallel if need be, but I want to know if there is a more efficient way to evaluate and plot a ParametricFunction over a range of values of the parameter?

It is faster to use NDSolve to obtain an interpolation over the parameter than solving the ODE over and over with different boundary conditions:

G = 0.00589667;
B = 2066.84;

ode = {Derivative[1, 0][y1][z, t] == 0.0039311133 y2[z, t] y3[z, t],
Derivative[1, 0][y2][z, t] == -2066.84 y3[z, t] - 0.0039311133 y3[z, t] Derivative[1, 0][y1][z, t],
Derivative[1, 0][y3][z, t] == 2066.84 y2[z, t]};

iv = {y1[0, t] == 0., y2[0, t] == 0., y3[0, t] == -Exp[-(t - 5)^2]};

sol = NDSolve[{ode, iv}, {y1, y2, y3}, {z, 0., 50.}, {t, 0., 10}];

Plot[Evaluate[y1[50, t] /. sol], {t, 0, 10}, PlotRange -> All]