Finite difference techniques can be applied to the numerical solution of the initial-boundary value problem in S for the semilinear Sobolev or pseudo-parabolic equation \[ \sum\limits_{i = 1}^n {\left[ {\frac{\partial } {{\partial x_i }}\left( {a_i \frac{\partial } {{\partial x_i }}u_t } \right) + \frac{\partial } {{\partial x_i }}\left( {b_i \frac{\partial } {{\partial x_i }}u} \right)} \right]} - q = ru_t , \] where $a_i $, $b_i $, q and r are functions of space and time variables, q is a boundedly differentiable function of u, and S is an open, connected domain in $\mathbb{R}^n $. Under suitable smoothness conditions, the solution of a Crank–Nicolson type of difference equation is shown to converge to u in the discrete $L^2 $-norm with an $O((\Delta x)^2 + (\Delta t)^2 )$ discretization error.The numerical problem is reduced to the inversion of a certain matrix at each time level. For the problem with constant coefficients in a two- or three-dimensional cube, a two-level iteration scheme with a Picard-...