I try to solve numerically the differential equation with mixed derivative. I found one possibility for the mixed derivative like y(i+1,j-1)-y(i,j-1)-y(i-1,j+1)+y(i-1,j). But may be there are another ways...
The discretisation you provided is incomplete. To compute the mixed derivative of u with respect to x and y you start from u_xy = (u_y)_x (or the other way around, your choice). Then you apply central differencing twice. The following applies to a uniform grid with spacing h. First find u_y at (i+1,j) and (i-1,j):
This discretisation is second order accurate. You will have to deal with the boundary points of your domain and may need to switch from central differences to forward or backward differences, depending on the boundary.
Also note that if u is a known function (rather than discrete values from an array for example), you can achieve much higher accuracy with Lagrange or Chebyshev polynomial approximations or using Ridders method of polynomial extrapolations.
Oops! I forgot to change the sign of the last term in the numerator. Copying and pasting is a bad practice.
R.C. Mittal's answer is for meshes that are uniform in x and y with hx and hy spacing in the x and y directions, respectively, hence the hx and hy. If you have a non-uniform mesh use the differences in the absissae and ordinates of the grid points as needed instead of h, hx or hy.
See my paper "Supraconvergence and supercloseness of a scheme for elliptic equations on non-uniform grids" Numerical Functional Analysis and Optimization, 27 (5-6). 539-564, 2006, where we propose a finite difference discretization using non-uniform grids and solutions with lower regularity smoothness.