We present a new methodology to perform the ǫ-expansion of hypergeometric functions with linear ǫ dependent Pochhammer parameters in any number of variables. Our approach allows one to perform Taylor as well as Laurent series expansion of multivariable hypergeometric functions. Each of the coefficients of ǫ in the series expansion is expressed as a linear combination of multivariable hypergeometric functions with the same domain of convergence as that of the original hypergeometric function thereby providing a closed system of expressions. We present illustrative examples of hypergeometric functions in one, two and three variables which are typical of Feynman integral calculus.