What is the maximum number of squares we can compact within a defined area? What is the maximum number of squares (with sides equal to $a$) that we can compact within a region limited by curves or lines defined by functions?
Note: So as not to make the problem more complex than it already is, consider squares perfectly coupled or aligned!
Let's see some examples:
The figure $1$ shows a model of how to arrange squares (side a) in rows and columns aligned in a region inside a parabola $f(x)$. Notice that we get only $12$ squares:

However, if we make a slightly more compact arrangement we can get a maximum of $13$ squares inserted! See figure $2$.
Furthermore the problem becomes more complex when we work with regions defined by two functions $f(x)$ and $g(x)$ as shown in figure $3$.
The question: How to determine the best arrangement of the squares and consequently the largest possible number of squares, using the maximum of the available area?
See also the link for a similar question asking the maximum area of a single square inserted in a defined region:
What is the maximum area of a square inserted within a region limited by certain functions?
 A: This is not the definitive answer to the most general problem, but it is a method for simpler cases, where we will make the following considerations:
The squares will be perfectly coupled and aligned, where rotations between them are not allowed.
For a better visualization and simplification, we analyze a convex region located in the $1st$ quadrant, delimited between $f(x)$, the $x$ axis and $x_{i} \leqslant x \leqslant x_{f}$, as shown in figure below.
Let us still consider simpler functions (curves), where they have $1$ or $2$ roots when we do $f(x) = na$, with $n \in Z_{+}^*$, and $a=$ side of square, where $a<y_{max}$ (maximum of $f(x)$ between $x_{i}$ and $x_{f}$).
See figure 1:

Let us consider simpler functions such as lines, parabolas, circumferences, etc, invertible [exist $f^{-1}(x)$] in defined interval, where we can more easily find the roots of the equations, when we equal $f(x)$ with the lines $na$.
$y = (na)$ are all parallel lines, where $n \in$ {$1, 2, 3, ..., n_{max}$}, with nmax given by:
$$n_{max} = \lfloor \frac{y_{max}}{a} \rfloor$$
For each range between $(na)$ and $(n-1)a$, delimited by the possible roots of the equation $f(x) = na$, we have the number of squares inserted in this range:
$$N_{na} = \lfloor \frac{\Delta_{na}}{a} \rfloor$$
Where $\Delta_{na} = X_{na}'' - X_{na}'$, with $X_{na}''$ and $X_{na}'$ are the roots of $ f(x) = na$ or yet $X_{na} = f^{-1}(na)$.
Thus, the total of squares inserted in the region is given by the sum of the maximum number of squares arranged in each range, that is:
$$N_{total} = \sum_{n=1}^{n_{max}} \lfloor \frac{\Delta_{na}}{a} \rfloor$$
We must also consider where the roots of the equations $f(x) = na$ are located with respect to the interval $H =$ {$x \in R |$ $x_{i} \leqslant x \leqslant x_{f}$}, and $S$ is the solution set of the equation:


*

*If $S =$ {$\emptyset$}, then $\Delta_{na} = X_{f} - X_{i}$

*If $S =$ {$x'$}, then:
$\Delta_{na} = X_{f} - X_{i}$, if $x'\leqslant X_{i}$ or $x'\geqslant X_{f}$
$\Delta_{na} = X_{f} - x'$, if $f(x)$ is growing in $x'$ and $x'\leqslant X_{f}$
$\Delta_{na} = x' - X_{i}$, if $f(x)$ is decreasing in $x'$ and $x'\geqslant X_{i}$


*

*If $S =$ {$x', x''$}, with $x''> x'$ then:


$\Delta_{na} = X_{f} - X_{i}$, if $x', x''\notin H$
$\Delta_{na} = X_{f} - x'$, if $x''\notin H$ and $x' \in H$
$\Delta_{na} = x'' - X_{i}$, if $x'\notin H$ and $x'' \in H$
$\Delta_{na} = x'' - X'$, if $x',x'' \in H$
Concave regions should be analyzed with greater care, as part of the range may be outside the region $R$, see point $P$ in figure below. Here $\Delta_{na}$ will be separated into two sub-regions:

In this case, for example: $\Delta_{na}' = x' - X_{i}$ and $\Delta_{na}'' = X_{f} - x''$
When we have a region between two functions (or curves) we can make the difference function between them $d(x) = f(x) - g(x)$ and apply the method to this resulting function, in the same interval? I'm not sure this can work for any functions, although the area of the region of the two graphs will be the same, see figure below:

Enough talk, let's go to an example: Let us determine the maximum number of squares of side $a = 1$, inserted in the inner region with a circle of radius = $5$, center at $(0,0)$, limited to the $1st$ quadrant:

By the feature of the function, we can easily see that all the roots of the $f(x) = na$ equation will be located in interval: $X_{i} \leqslant x \leqslant X_{f}$, with  $X_{i}=0$ and $X_{f}=5$.
We also have to:
$a = 1$
$n_{max} = 5$
$f^{-1}(x) = \sqrt{25-x^2}$
$X_{na} = f^{-1}(na)$
$\Delta_{na} = X_{na} - X_{i}$, so $\Delta_{n} = X_{na} - 0 = X_{n}= f^{-1}(n)$
Therefore:
$$N_{total} = \sum_{n=1}^{5} \lfloor \frac{\Delta_{n}}{1} \rfloor = \sum_{n=1}^{5} \lfloor f^{-1}(n) \rfloor = \sum_{n=1}^{5} \lfloor \sqrt{25-n^2} \rfloor$$
So:
$$N_{total} = \lfloor \sqrt{25-1^2} \rfloor + \lfloor \sqrt{25-2^2} \rfloor + \lfloor \sqrt{25-3^2} \rfloor + \lfloor \sqrt{25-4^2} \rfloor + \lfloor \sqrt{25-5^2} \rfloor$$
$$N_{total} = 4 + 4 + 4 + 3 + 0$$
$$N_{total} = 15$$
I think this may help to analyze a more general question.
