First define a symbolic function f(u1,u2,...,ui) using the syms command.
Next, define a vector u=[u1,u2,...,ui] containing all variables of f.
Define an 1 x i row vector q=[q1,q2,...,qi], containing the values that the function will be evaluated in. Define also an sq=[sq1,sq2,...,sqi] vector, containing the uncertainty on every value in q.
q and sq can also be
m x i matrices (for varying values of qi and sqi). You can also have varying qi and a fixed sqi (the function will recognize if sqi is fixed (a row vector) or changing (a matrix).
The function will output an m x 2 matrix with the first column containing f(q) and the second column the error propagation of f(q) using this formula: