  The hypergeometric function method naturally provides the analytic
expressions of scalar integrals from concerned Feynman diagrams in some
connected regions of independent kinematic variables, also presents the systems
of homogeneous linear partial differential equations satisfied by the
corresponding scalar integrals. Taking examples of the one-loop $B_{_0}$ and
massless $C_{_0}$ functions, as well as the scalar integrals of two-loop vacuum
and sunset diagrams, we verify our expressions coinciding with the well-known
results of literatures. Based on the multiple hypergeometric functions of
independent kinematic variables, the systems of homogeneous linear partial
differential equations satisfied by the mentioned scalar integrals are
established. Using the calculus of variations, one recognizes the system of
linear partial differential equations as stationary conditions of a functional
under some given restrictions, which is the cornerstone to perform the
continuation of the scalar integrals to whole kinematic domains numerically
with the finite element methods. In principle this method can be used to
evaluate the scalar integrals of any Feynman diagrams.
