A new global spatial discretization method is developed to accurately calculate natural frequencies and dynamic responses of two-dimensional continuous systems such as membranes and Kirchhoff plates. The transverse displacement of a two-dimensional continuous system is separated into a two-dimensional internal term and a two-dimensional boundary-induced term; the latter is interpolated from one-dimensional boundary functions that are further divided into one-dimensional internal terms and one-dimensional boundary-induced terms. The two- and one-dimensional internal terms are chosen to satisfy predetermined boundary conditions, and the two- and one-dimensional boundary-induced terms use additional degrees of freedom at boundaries to ensure satisfaction of all boundary conditions. A general formulation of the method that can achieve uniform convergence is established for a two-dimensional continuous system with an arbitrary domain shape and arbitrary boundary conditions, and it is elaborated in detail for a general rectangular Kirchhoff plate. An example of a rectangular Kirchhoff plate that has three simply-supported boundaries and one free boundary with an attached Euler-Bernoulli beam is investigated using the developed method and results are compared with those from other global and local spatial discretization methods. Natural frequencies and dynamic responses that include the displacement, the velocity, rotational angles, a bending moment, and a transverse shearing force are calculated using both the developed method and the assumed modes method, and compared with results from the finite element method and the finite difference method, respectively. Advantages of the new method over local spatial discretization methods are much fewer degrees of freedom and much less computational effort, and those over the assumed modes method are better numerical property, a faster calculation speed, and much higher accuracy in calculation of the bending moment and the transverse shearing force that are related to high-order spatial derivatives of the displacement of the plate with an edge beam.