After seeing Ryan Budney's function, I came up with this:
On an n-torus $\mathbb{T}^n$ we have a Morse function $f=\sum^n_{i=1}\cos\theta_i$, which has critical points $\text{crit}(f)=\lbrace(\theta_1,\ldots,\theta_n)\;|\;\theta_i = 0\text{ or }\pi\;\;\forall\;i\rbrace$.
In particular, we get a self-indexing Morse function $f(\theta_1,\theta_2,\theta_3)=\frac{1}{2}\sum^3_{i=1}\cos\theta_i+\frac{3}{2}$ on $\mathbb{T}^3$.
Indeed, it has a unique maximum $(0,0,0)$ with value 3, and a unique minimum $(\pi,\pi,\pi)$ with value 0. There are 3 index-1 critical points (two angles are $\pi$ and one angle is $0$) with value $1$, and 3 index-2 critical points (one angle is $\pi$ and two angles are $0$) with value $2$, and the indices are evident from the Hessian $H_{(\theta_1,\theta_2,\theta_3)}f=-\text{diag}(\cos\theta_1,\cos\theta_2,\cos\theta_3)$.
Now this $f$ then provides a Heegaard diagram for $\mathbb{T}^3$ of genus $3$ (since there are 3 index-1 and 3 index-2 critical points, corresponding to the $\alpha,\beta$-curves). As Wikipedia states: It was proved by Frohman and Hass that any other genus 3 Heegaard splitting of the three-torus is topologically equivalent to this one (the one in the original question).
Is it true then that my $f$ and Ryan's $f$ provide an example of the "gap" $Diff(\mathbb{T}^3,\mathcal{H})$ between Heegaard decompositions and self-indexing Morse functions, as explained in Ryan's comment?