This paper introduces multidimensional block-pulse functions and develops a method of numerically integrating a class of partial differential equations based on orthogonal approximation of a function of several variables. The resulting solutions are piece-wise constant with minimal mean square error. The original partial differential equation is transformed into a computationally convenient algebraic form. Tha new operational matrices for partial operators presented in this paper are structurally related to the corresponding operational matrices of ordinary-calculus (i.e. the calculus of a single variable) and appear as Kronecker products. Certain algebraic properties of these matrices suggest enormous reduction in the computational effort associated with the algebra of these apparently large and sparse matrices. All solutions are reduced to simple recurrence form and the stability of these recurrence solutions is investigated. Examples are included to illustrate the proposed method in the study of distributed parameter systems.