I've been trying to answer this question for several years and it turned out to be really hard, even for the $2$-sphere. Below I will discuss this case.
First of all one should ask what is the number $m(k)$ of topological types of (stable) Morse functions on $S^2$ with precisely $k$ saddle points. (such a function has $2k+2$ critical points.) I showed that the generating series
$$ x(t) := \sum_{k\geq 0} \frac{m(k)}{(2k+1)!} t^{2k+1}, $$
is the inverse of an elliptic integral; see this paper. More precisely $x(t)$ is the inverse of the function
$$ x\mapsto t(x)=\int_0^x \frac{ds}{\sqrt{s^4/4-s^2-2sx+1}} ds. $$
This fact leads to a positive answer to a question of V.I. Arnold who conjectured that
$$\log m(k)\sim 2k\log k $$
as $k\to \infty $. I refer you to this paper for details. This shows that $m(k)$ grows rather fast as $k\to \infty$.
Any polynomial $P$ of degree $d$ in $\newcommand{\bR}{\mathbb{R}}$ on $\bR^n$ can be uniquely decomposed as a sum
$$ P= \sum_{0\leq j+2k\leq d} r^{2k} H_{j}, \;\; r^2= (x_1^2+\cdots +x_n^2), $$
where $H_{j}$ is a darmonic polynomial of degree $j$. On $\bR^3$ the space of degree $d$ hormonic polynomials has dimension $2d+1$. If we denote by $U_d$ the subspace of $C^\infty(S^2)$ consisting of the restrictions to $S^2$ of the polynomials of degree $\leq d$ we deduce that
$$\dim U_d=\sum_{0\leq k\leq d} (2k+1)=(d+1)^2. $$
Denote by $K_d$ the expected number of critical points of a random function in $U_d$. I showed that
$$ K_d\sim C\dim U_d\sim Cd^2 $$
as $d\to \infty$ where $C$ is a certain explicit constant; see this paper and this paper.
It turns out that the number of critical points of a random function in $U_d$ is highly concentrated around its mean $K_d$, i.e., the probability that the number of critical points of a random function in $U_d$ is far from the mean $K_d$ is extremely small as $d\to\infty$. In more precise technical terms, the variance of the (random) number of critical points of a (random) function in $U_d$ has the same size as $K_d$, which makes the standard deviation of size $\sqrt{K_d}$, much, much smaller than $K_d$ for $d$ large.
I personally believe, based on some empirical evidence, that the mean is close to the maximum number of critical points in the sense that if we denote by $\mu_d$ the maximum number of critical points of a Morse function in $U_d$, then $\mu_d \sim C'' d^2$ as $d\to\infty$.
My guess is that the number of topological types of functions in $U_d$ as $d\to \infty$ is roughly
$$ \sum_{k=1}^{K_d/2} m(k), $$
where I recall that $m(k)$ denotes the number of topological types of Morse functions with $k$ saddle points, i.e., $2k+2$ critical points.
