Your question is a bit too vague. Or rather, the input-data you have is a bit vaguely described. The answer depends on how you'll specify that input...
If you have a reasonably low number of particles and their trajectories (something like 
) you could simply go for a looped version of plot3, maybe something like:
) you could simply go for a looped version of plot3, maybe something like:clf
hold on
E_kin = squeeze(sum(v(:,1,:).^2,3));
g = (E_kin-min(E_kin))/(max(E_kin)-min(E_kin));
r = (r(:,1,1)-min(r(:,1,1)))/(max(r(:,1,1))-min(r(:,1,1)));
b = (r(:,1,2)-min(r(:,1,2)))/(max(r(:,1,2))-min(r(:,1,2)));
for i1 = num_particles:-1:1
ph(i1) = plot3(r(i1,:,1),r(i1,:,2),r(i1,:,3));
set(ph(i1),'color',[r(i1),g(i1),b(i1)])
end
That would plot the trajectories of all your particles in 3-D with the trajectories coloured according to their
initial x-y position (red and blue), and initial kinetic energy (green).
On the other hand that seems like such an obvious solution, so you might such a large number of particles that this doesnt work. For that perhaps you can use local averaged velocities, but that would involve a lot of looping:
for each time-step for each sub-volume, find all all particles inside, calculate their average velocity, when done with all subvolumes, you should be able to use streamline3, then on to next time-step
HTH
