  Starting from the Mellin-Barnes integral representation of a Feynman integral
depending on set of kinematic variables $z_i$, we derive a system of partial
differential equations w.r.t.\ new variables $x_j$, which parameterize the
differentiable constraints $z_i=y_i(x_j)$. In our algorithm, the powers of
propagators can be considered as arbitrary parameters. Our algorithm can also
be used for the reduction of multiple hypergeometric sums to sums of lower
dimension, finding special values and reduction equations of hypergeometric
functions in a singular locus of continuous variables, or finding systems of
partial differential equations for master integrals with arbitrary powers of
propagators. As an illustration, we produce a differential equation of fourth
order in one variable for the one-loop two-point Feynman diagram with two
different masses and arbitrary propagator powers.
