Subscript the variables so they are column vectors, then put them in the table and return the table as one of the outputs (or perhaps the only output).
dx(:,it) = eliminacao_gauss( J(xo), -f(xo) );
erro(:,it) = norm(xk-xo)/max(abs(xk));
s = table(itv,dx,xk,erro);
Note that ‘s’ is not otherwise defined anywhere, so I assigned the table array to it.