First, a note about those starting position plots:
"Each column, 1-3, corresponds to one of the three bodies and the rows is that x,y, or z positions as time changes."
Then the starting position of the first body is
[x(1,1) y(1,1) z(1,1)]
(You have this correct in your plot.)
But the starting position of the second body is
[x(1,2) y(1,2) z(1,2)]
not
[x(2,2) y(2,2) z(2,2)]
which is what you have in your plot, and which is actually the position of the second body at time 2. (Similar for the third body / time 3.)
Now, to answer your question about the center of mass:
center_of_mass = [mean(x,2) mean(y,2) mean(z,2)]
mean(_,2) takes the mean of each row (dimension 2) of the given array. In this case, that's the mean of each body's x, then y, then z (for all rows of each matrix), concatenated together in a single matrix called center_of_mass. (Or you could define 3 separate variables for center of mass x-coordinate, center of mass y-coordinate, center of mass z-coordinate, but I won't do that.)
Then to plot the center of mass:
plot3(center_of_mass(:,1),center_of_mass(:,2),center_of_mass(:,3))
(Or if you just want to calculate it and plot it without storing it:)
plot3(mean(x,2),mean(y,2),mean(z,2))
