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
$$ f(t) := \sum_{k\geq 0} \frac{m(k)}{(2k+1)!} t^{2k+1}, $$
is the inverse of an elliptic integral; see this paper. 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. $$
I showed that as $d\to \infty$ the expected number of critical points of a function in $U_d$ is comparable to $C \dim U_d\sim C'd^2$, where $C, C'$ are certain explicit constants; see this paper and this paper. Denote by $K_d$ this expected number. Moreover 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 a random function in $U_d$ a number of critical far from the mean $K_d$ is extremely small as $d\to\infty$.
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 $m(k)$ I recall that $m(k)$ denotes the number of topological types of Morse functions with $k$ saddle points, i.e., $2k+2$ critical points.