Simple grid over geographic coordinates

I wanted to build a simple grid over the earth’s surface. I need this functionality for a computer program where I quickly need to check wheter a given coordinate in (latitude, longitude) is inside a cluster.
These cluster are therefore represented as index-pairs (latitude-index, longitude-index). Which I was computing quite naively, independent in each dimension. Here is some Java code on how I compute these indices:

public GeoIndex getIndex(double latitude, double longitude, double cellSize) {
    if (!isPointInBoundaries(latitude, longitude)) {
        return null;
    }
    long latitudeIndex = (long) Math.floor(haversine(boundaries.latBottomLeft, boundaries.longBottomLeft,
                                                    latitude, boundaries.longBottomLeft) / cellSize);
    long longitudeIndex = (long) Math.floor(haversine(boundaries.latBottomLeft, boundaries.longBottomLeft,
                                                    boundaries.latBottomLeft, longitude) / cellSize);

return new GeoIndex(latitudeIndex, longitudeIndex);
}

Here the cellSize is supposed to be the grid length for each grid-cell. And the boundaries object gives a point that is the 0-point in my coordinate system, so to say. The lowest, left-most point from where the coordinate system will be originated. I want to cover regions like, Europe, North America, China…

First question, is this a solid approach and will points always map to the same grid cells? Or is this too simplistic.

My real issue started when I was trying to compute the boundary points of each grid-cell. I though I could just recompute the boundaries. But it seems to be an issue that the earth radius shrinks the farther north you go. Is there a way to recompute the boundaries of each cell? I was trying the following approach:

static toCoordinates(geoIndex) {
    const coordinates = [];
    const cells = geoCluster.getAiActivatedMapCells();
    const distanceInKm = cellSize / 1000.0;
    const latitudePoint = GeoUtils.calculateEndPointWithDistance(
        boundaries.latBottomLeft,
        boundaries.longBottomLeft,
        distanceInKm * geoIndex.latitude,
        0.0
      );
      const longitudePoint = GeoUtils.calculateEndPointWithDistance(
        boundaries.latBottomLeft,
        boundaries.longBottomLeft,
        distanceInKm * geoIndex.longitude,
        90.0
      );
      const bottomLeftPoint = {
        lat: latitudePoint.lat,
        lng: longitudePoint.lng,
      };

      const topLeftPoint = GeoUtils.calculateEndPointWithDistance(
        latitudePoint.lat,
        longitudePoint.lng,
        distanceInKm,
        0.0
      );
      const bottomRightPoint = GeoUtils.calculateEndPointWithDistance(
        latitudePoint.lat,
        longitudePoint.lng,
        distanceInKm,
        90.0
      );
      const topRightPoint = {
        lat: topLeftPoint.lat,
        lng: bottomRightPoint.lng,
      };

      return {
        topLeft: topLeftPoint,
        topRight: topRightPoint,
        bottomRight: bottomRightPoint,
        bottomLeft: bottomLeftPoint,
      };
}

static calculateEndPointWithDistance(
    startLatitude,
    startLongitude,
    distanceInKm,
    bearing
  ) {
    const distanceRatio = distanceInKm / GeoUtils.EARTH_RADIUS;
    const startLat = GeoUtils.deg2rad(startLatitude);
    const startLon = GeoUtils.deg2rad(startLongitude);

    const endLat = Math.asin(
      Math.sin(startLat) * Math.cos(distanceRatio) +
        Math.cos(startLat) *
          Math.sin(distanceRatio) *
          Math.cos(GeoUtils.deg2rad(bearing))
    );
    const endLon =
      startLon +
      Math.atan2(
        Math.sin(GeoUtils.deg2rad(bearing)) *
          Math.sin(distanceRatio) *
          Math.cos(startLat),
        Math.cos(distanceRatio) - Math.sin(startLat) * Math.sin(endLat)
      );

    return { lat: GeoUtils.rad2deg(endLat), lng: GeoUtils.rad2deg(endLon) };
}

Sorry that this is JavaScript. But basically I try to recompute the latitude and longitude of the lower left point of the grid cell independently again. With the function calculateEndPointWithDistance which is an implementation of the the function I found here:
https://www.movable-type.co.uk/scripts/latlong.html (Destination point given distance and bearing from start point). And the compute the grid boarders.

However, the longitude is always way off. I don’t really understand way. I understand why the grid cells are to big and overlap, but I don’t understand why they are way off to west at all. Can you explain me why?