The paper presents a step-by-step construction of a finite-difference scheme for a heterogeneous biharmonic equation under zero boundary conditions superimposed on the desired function and its first-order partial derivatives. The finite-difference scheme is based on a square twenty-five-point pattern and has an implicit character. On analytical grid, the error of approximation of the biharmonic operator by the difference analog and the error of approximation of boundary conditions imposed on the first order partial derivatives are calculated by the expansion of the function in the Taylor series with the remainder term in the form of a Lagrange. The boundary conditions imposed on the sought function are satisfied precisely. A finite-difference scheme approximates a boundary value problem with a second order of accuracy over the mesh step. With the help of the Maple computer algebra system the solutions of the problem for different grid steps are obtained. The dependence of the minimum function and calculation time on the number of significant digits is revealed. The optimal number of significant digits is found. The convergence rate of the numerical scheme is analyzed. The dependence of the minimum value of the function and the calculation time on the value of the grid step is established. The optimal step value is found. A three-dimensional graph of the solution and its profiles in the middle sections are constructed. The advantages of the developed finite-difference scheme are indicated. Obtained results correspond to the physical meaning of the problem and are consistent with similar numerical and approximate analytical solutions.