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][1].  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][2]   who conjectured that  
$$\log m(k)\sim 2k\log k $$

 as $k\to \infty $. I refer you to  [this paper][3] 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][4] and [this paper][5].   

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.


  [1]: http://www3.nd.edu/~lnicolae/Morse-count.pdf
  [2]: http://www.mathnet.ru/php/presentation.phtml?option_lang=eng&presentid=189
  [3]: http://www3.nd.edu/~lnicolae/statistics.pdf
  [4]: https://arxiv.org/abs/1008.5085
  [5]: http://www3.nd.edu/~lnicolae/CritSetStat.pdf