I suggest to rely on a recursive solution like that in the attached spreadsheet. It is based on the Vincenty's formula which provides accurate avaluation (1/10^-6) of the distance between two points of given lat and long (dec degrees).

instructions:

Allow macro to run.

provide lat and lon in cells a3 and b3 respectively.

provide the step (dec degrees)

Push the green button.

read area (m^2) with the desired approximation .

there are many sites in which you can get the Vincenty's algorithm to adapt it to your needs. My solution is but an example.

Dear Pietro,

Thank you very much for your help. Using your indications, I have found an script for Matlab very useful.

Thanks again.

http://www.mathworks.es/es/help/map/ref/areaint.html

Very simple: it's just the angular part of the volume element in spherical polar coordinates exchanging latitude for polar angle: if a is longitude and b is latitude, dOmega = cos(b) db da.

You can also apply the Gauss-Bonnet formula to the "spherical rectangle". The "parallels"(latitudes) are not geodesics, but they have constant geodesic curvature (a simple trigonometric function depending on the latitude)