# Trouble with ParametricNDSolveValue

I have this:

eqns = {
x1'[t] == -10 x1[t] + 10 x2[t],
x2'[t] == a1 x1[t] - x2[t] - x1[t]*x3[t],
x3'[t] == -(8/3) x3[t] + x1[t] x2[t]
};
ics = {x1 == a, x2 == b, x3 == c};


Then this:

pfun = ParametricNDSolveValue[{eqns, ics}, {x1, x2, x3}, {t, 0,
10}, {a1, a, b, c}]


But when I try to plot:

Plot[pfun[50, 1, 1, 1][t], {t, 0, 10}]


I get a blank image. What am I doing incorrectly?

Due to Helpful Answers: Thanks for the responses. Turns out the easiest answer for me is to replace the output {x1, x2, x3} with {x1[t], x2[t], x3[t]}.

pfun = ParametricNDSolveValue[{eqns, ics}, {x1[t], x2[t], x3[t]}, {t, 0,
10}, {a1, a, b, c}];
Plot[Evaluate[pfun[50, 1, 1, 1]], {t, 0, 10}]


Which produces this image: You can see why this works by entering:

pfun[50, 1, 1, 1]


Which produces this output:

{InterpolatingFunction[{{0., 10.}}, <>][t], InterpolatingFunction[{{0., 10.}}, <>][t], InterpolatingFunction[{{0., 10.}}, <>][t]}

I can also grab the first solution for x1 with pfun[50, 1, 1, 1][], which allows me to do things like the following:

Plot[Evaluate[{pfun[50, 1, 1, 1][], pfun[50, 1.1, 1, 1][]}], {t,
0, 10}]


Which produces the following image: See some further amazing help here. Thanks, everyone.

Plot[Evaluate[Through@pfun[50, 1, 1, 1]@t], {t, 0, 10}]
(* or Plot[#[t] & /@ pfun[50, 1, 1, 1], {t, 0, 10}, Evaluated -> True] *) Note: The parametric function pfun[50, 1, 1, 1] is a list of 3 InterpolatingFunctions:

pfun[50, 1, 1, 1]


{InterpolatingFunction[{{0., 10.}}, <>], InterpolatingFunction[{{0., 10.}}, <>], InterpolatingFunction[{{0., 10.}}, <>]}

To evaluate these functions at x you need to use

Through@pfun[50, 1, 1, 1]@x


or

#[x] & /@ pfun[50, 1, 1, 1]


Both return

{InterpolatingFunction[{{0., 10.}}, <>][x], InterpolatingFunction[{{0., 10.}}, <>][x], InterpolatingFunction[{{0., 10.}}, <>][x]}

which can be wrapped with Evaluate and used as the first argument of Plot. Alternatively, you can use the Plot option Evaluated->True without the Evaluate wrapper.

The combination,

pfun = ParametricNDSolveValue[{eqns, ics}, {x1[t], x2[t], x3[t]}, {t, 0, 10},
{a1, a, b, c}];
Plot[Evaluate[pfun[50, 1, 1, 1]], {t, 0, 10}]


also works, giving the same curves.