Animate flow lines of time-dependent 3D dynamical system

You can use ParametricPlot3D to get smoother orbits:

SeedRandom[1]
n = 50;
seeds = RandomReal[{-1, 1}, {n, 3}];
tbar = 10;

Animate[ParametricPlot3D[Evaluate[func[seeds[[#]], t] & /@ Range[Length@seeds]], 
   {t, 0, tmax}, 
    BoxRatios -> 1, 
    PlotStyle -> Arrowheads[Medium],
    ImageSize -> 400,
    PlotRange -> {{-60, 60}, {-1, 1}, {-10, 10}}] /. Line -> Arrow,
  {tmax, .01, tbar}]

enter image description here

Update: "Is there a way to use this this implementation to animate as all of orbit 1, followed by all of orbit 2, followed by....?"

colors = Table[Hue@RandomReal[], {n}];

Animate[ParametricPlot3D[Evaluate[ConditionalExpression[func[seeds[[#]], t - (#-1) tbar],
      (# - 1) tbar <= t <= # tbar] & /@ Range[n]], 
   {t, 0, tmax}, 
   BaseStyle -> Arrowheads[Medium], 
   BoxRatios -> 1, 
   ImageSize -> 400, 
   PlotStyle -> colors,
   PerformanceGoal -> "Quality", 
   PlotRange -> {{-50, 60}, {-1, 1}, {-5, 15}}] /. Line -> Arrow, 
 {tmax, .01, n tbar}, 
 AnimationRate -> 10]

enter image description here


Not exactly what was asked, but another way to visualize the flow, based on How can I create a fountain effect?:

DynamicModule[
 {x0, y0, z0, last = 0, lam = 1.5, n = 500, colors, replace},
 last = Clock[Infinity];
 {x0, y0, z0} = RandomReal[{-1, 1}, {3, n}];
 colors = RandomColor[n];

 Graphics3D[GraphicsComplex[
   Dynamic@ With[{t = Clock[Infinity]},
     With[{dt = (t - last)/2}, With[{dl = lam^dt},
       last = t;
       x0 = x0*dl; y0 = y0/dl; z0 = z0 + dt; (* integration of velocity *)
       replace = Pick[Range@n, UnitStep[z0 - 1], 1];
       x0[[replace]] = RandomReal[{-1, 1}, Length@replace];
       y0[[replace]] = RandomReal[{-1, 1}, Length@replace];
       z0[[replace]] = RandomReal[{-1, -1 + dt}, Length@replace];
       Transpose@{x0, y0, z0}
       ]]],
   Point[Range@n, VertexColors -> colors]], 
  PlotRange -> {{-2, 2}, {-2, 2}, {-1, 1}}, Axes -> True, 
  AxesLabel -> {x, y, z}],

 Initialization :> ({x0, y0, z0} = RandomReal[{-1, 1}, {3, n}])
 ]

enter image description here

The integration is based on the ODE for $\Phi$, which is autonomous and linear and can be done by scalings and translation: $${d \over dt}\,(x,y,z) = (x \log \lambda, -y \log \lambda, 1)$$