A Nyström method for second-kind Volterra integral equations on the square

In this paper, we propose a Nyström-type method for the numerical approximation of second-kind bivariate Volterra integral equations on the unit square. The method is based on bivariate Generalized Bernstein (GB) polynomials and on a cubature formula defined on equally spaced nodes, allowing the integral operator to be discretized directly without requiring any change of variables. The use of uniform grids makes the proposed approach particularly suitable for applications in which data are naturally sampled at equidistant points.The approximation framework relies on tensor-product GB polynomials whose degrees may differ with respect to each variable. This flexibility allows different approximation orders to be prescribed according to the regularity of the target function in each variable. As a result, the computational cost can be reduced while preserving high level of accuracy. Compared with classical Bernstein operators, the proposed GB operators achieve higher approximation orders depending on the smoothness of the approximated functions.The resulting Nyström scheme leads to a well-conditioned linear system and is proved to be stable and convergent. Error estimates are derived in Sobolev and Zygmund spaces. Several numerical experiments are presented to validate the theoretical results and to illustrate the effectiveness of the method for kernels and right-hand side functions with different smoothness properties.