As far as I can tell, there is no way to create a quiver plot (or quiver3 plot) with quivers of different colors.
One way to approach this problem is to create one quiver3 plot (of a given color) for each color in the colormap (i.e., the colors depicted in the colorbar). Now, in quiver/3 the quivers are scaled in length so as not to overlap each other (with an additional optional scale factor), but the only way (as far as I can tell) to get multiple quiver/3 plots to use the same scale factor is to turn that scaling off. So that could work, scaling the magnitude data down by an appropriate amount before sending it to quiver/3:
xyz = readmatrix('vector3dn.txt');
xyz(:,4:6) = xyz(:,4:6)*sf;
magnitude = sqrt(sum(xyz(:,4:6).^2,2));
floor(min(magnitude)/mag_res) ...
ceil(max(magnitude)/mag_res) ...
mthresholds = linspace(mlim(1),mlim(2),n_colors+1);
idx = magnitude >= mthresholds(ii) & magnitude < mthresholds(ii+1);
args = num2cell(xyz(idx,:),1);
quiver3(args{:},'off','Color',cmap(ii,:))
Or you can use line objects instead of quivers; the approach is the same - splitting it up into multiple sets of lines, one for each color. You just have to calculate the "other" end of each line (x+u, y+v, z+w).
xyz = readmatrix('vector3dn.txt');
magnitude = sqrt(sum(xyz(:,4:6).^2,2));
floor(min(magnitude)/mag_res) ...
ceil(max(magnitude)/mag_res) ...
mthresholds = linspace(mlim(1),mlim(2),n_colors+1);
idx = magnitude >= mthresholds(ii) & magnitude < mthresholds(ii+1);
x = xyz(idx,1)+[zeros(n,1) sf*xyz(idx,4)];
y = xyz(idx,2)+[zeros(n,1) sf*xyz(idx,5)];
z = xyz(idx,3)+[zeros(n,1) sf*xyz(idx,6)];
line(x.',y.',z.','Color',cmap(ii,:))