In this paper the one-and two-dimensional pseudoparabolic equations with nonlocal boundary conditions are approximated by the Euler finite difference scheme. In the case of classical boundary conditions the stability of all schemes is investigated by the spectral method. Stability regions of finite difference schemes approximating pseudoparabolic problem are compared with the stability regions of the classical discrete parabolic problem. These results are generalized for problems with nonlocal boundary conditions if a matrix of the finite difference scheme can be diagonalized. For the two-dimensional problem an efficient algorithm is constructed, which is based on the combination of the FFT method and the factorization algorithm. General stability results, known for the three level finite difference schemes, are applied to investigate the stability of some explicit approximations of the two-dimensional pseudoparabolic problem with classical boundary conditions. A connection between the energy method stability conditions and the spectrum Hurwitz stability criterion is shown. The obtained results can be applied for pseudoparabolic problems with nonlocal boundary conditions, if a matrix of the finite difference scheme can be diagonalized.