This paper presents an approach to the two-dimensional analysis of elastic isotropic deep beams using the finite difference method (FDM). Deep beams are subjected to in-plane loading and present a shear span to height ratio of less than 2.50; consequently, Euler-Bernoulli beam theory and Timoshenko beam theory do not apply. Deep beams analysis is generally conducted using numerical methods such as the finite element method and to a lesser extent the FDM; the strut-and-tie model and the stress field method are also widely utilized. Analytical approaches usually make use of the Airy stress function, where stresses are formulated in terms of the stress function; however, the exact solution of this function satisfying all of the boundary conditions can hardly be found, even for simple cases. In this paper, deep beams were analyzed using the FDM. The FDM is an approximate method for solving problems described with differential equations. The FDM does not involve solving differential equations; equations are formulated with values at selected nodes of the structure. Therefore, the deep beam was discretized with a two-dimensional grid, and additional nodes were introduced at the boundaries and at positions of discontinuity (openings, brutal change of material properties, non-uniform grid spacing), the number of additional nodes corresponding to the number of boundary conditions at the node of interest. The introduction of additional nodes allowed us to apply the governing equations at boundary nodes and satisfy the boundary and continuity conditions. An Airy stress function approach and a displacement potential function approach were considered in this study whereby strong formulations of equations (equilibrium, kinematic, and constitutive) were set. Stress and stability analyses were carried out with this model; furthermore, deep beams of varying stiffness, layered beams, and beams having openings were analyzed. For slender beams, the results obtained with the Airy stress function approach showed good agreement with those of the Euler-Bernoulli beam theory, and for deep beams, the shapes of stress distributions were in good agreement with a proper understanding of the behavior of structures. On the other hand, the displacement potential function approach delivered unsatisfactory results, probably due to the use of an inefficient equation solver; a more powerful tool will be needed in future research for this purpose.