def sum_squared_xy_derivative(xy_point, coeff_mat):
    """For a coordinate, and function, finds df/dx and df/dy and returns the
    sum of the squares

    Arguments:
        xy_point (tuple): (x,y)
        coeff_mat (np.array): Matrix of coefficients of the n order polynomial

    Returns:
        (float): (df/dx + df/dy)^2 where at a stationary point ~ 0
    """
    order = coeff_mat.shape[0]
    x, y = xy_point
    dx, dy = 0, 0

    for i in range(order):  # x index
        for j in range(order):  # y index
            if i > 0:
                dx += coeff_mat[i, j] * i * x**(i-1) * y**j
            if j > 0:
                dy += coeff_mat[i, j] * x**i * j * y**(j-1)

    return dx**2 + dy**2