The question has a negative answer if $M$ is compact. The reason for this is that any smooth function on a compact manifold has at least one critical point.
However, if you impose a weaker nondegeneracy condition, namely that the above matrix is nondegenerate only along the diagonal of $M\times M$, then the question has a positive answer and some interesting geometric consequences. Here are the details.
For any $f\in C^\infty(M\times M$ denote by $B_f(x,y)$ the matrix you wrote viewed as a bilinear map $\newcommand{\bR}{\mathbb{R}}$
$$ T_xM\times T_yM\to \bR. $$
Equivalently we can view this as a linear map
$$T_xM\to T^*_y M . $$
Note that
$$B_{f+g}= B_f+ B_g. $$
Seek $f$ as a sum of functions of the from $u(x)v(y)$, $u,v\in C^\infty(M)$. The vector space $\newcommand{\eF}{\mathscr{F}}$ $\eF$ of functions of this type is dense in $C^\infty(M\times M)$. Assume from now on that $M$ is compact. Set $m:=\dim M$.
If there existed a function $f\in C^\infty(M)$ such that $B_f(x,y)$ is nondegenerate for any $x,y\in M$, then you could find $h\in \eF$ such that $B_h(x,y)$ is nondegenerate for any $x,y\in M$.
Observe that
$$ B_{uv}(x,y)= du(x)\otimes dv(y).$$
The linear map $B_{uv}:T_xM\to T^*_y M$ has the form
$$ T_xM\ni X\mapsto (Xu)(x) dv(y)\in T_y^* M. $$
Its adjoint $B_{uv}^*: T_yM\to T_x^*M$ is also nondegenerate and it is given by
$$ T_yM\ni Y\mapsto (Yv)(y)du(x)\in T_x^* M. $$
Suppose that
$$f=\sum_{k=1}^N u_k(x)v_k(y). $$
Then
$$ B_f=\sum_{k=1}^N du_k(x)\otimes dv_k(y) $$
If $B_f(x,y)$ is nondegenerate, we deduce that $B_f(x,y): T_xX\to T^*_y M$ is onto and we conclude that the the covectors $dv_1(y),\dotsc, dv_N(y)$ span $T_y^* M$, for any $y$. Similarly, the covectors $du_1(x),\dotsc, du_N(x)$ span $T_xM$ for any $x\in M$. These two conditions can be conveniently rephrased as follows: the smooth maps $ \Phi,\Psi: M\to\bR^N$ defined by
$$\Phi(x)= \bigl(\; u_1(x),\dotsc, u_N(x)\;\bigr),\;\; x\in M, $$
$$\Psi(y)=\bigl(\; v_1(y),\dotsc, v_N(y)\;\bigr),\;\;y\in M $$
are immersions.
Fix a finite dimensional vector Euclidean space $E$ (think of $E$ as the above $\bR^N$) and immersions
$$\Phi, \Psi: M\to E. $$
For every $x\in M$ we have inclusions
$$i_{\Phi,x}, i_{\Psi,x}: T_x M\to E, $$
and surjections
$$ P_{\Phi,x}, P_{\Psi,x}: E^*\to T_x^* M. $$
More precisely $i_{\Phi_x}= D_x\Phi$, $P_{\Phi_,x}= (D_x\Phi)^*$, where $D_x\Phi$ is the differential of $\Phi$ at $x\in M$. If $\DeclareMathOperator{\Gr}{\boldsymbol{Gr}}$ if $\Gr_m(M)$ denotes the Grassmannian of $m$-dimensional subspaces of $E$, then the maps
$$ M\ni X \mapsto D_x\Phi(T_xM), D_x\Psi(T_xM)\in \Gr_m(E) $$
are the Gauss maps of the immersions $\Phi$ and $\Psi$.
Using the metric on $E$ we can regard $P_{\Phi,x}$ as a surjection $E\to T_xM$, the orthogonal projection onto the subspace $T_xM\subset E$.
In the compact case the problem you ask is now reduced to the following problem.
$\mathbf{P_1}$ Find immersions $\Phi,\Psi: M\to E$, $E$ Euclidean space such that for any $x,y \in M$ the operator $$ P_{\Psi_y}\circ i_{\Phi_x}=(D_y\Psi)^*D_x\Phi: T_xM \to T_yM $$ is onto?
Denote by $g: E\times E\to\bR$ denotes the inner product on $E$. Problem $\mathbf{P_1}$ is equivalent to the following.
$\mathbf{P_2}$ Find immersions $\Phi,\Psi: M\to E$ such that, if $$ f:=(\Phi\times \Psi)^* g, $$ then $B_f(x,y)$ is nondegenerate for any $x,y\in M$.
The last formulation shows that the example $f=x\cdot y$ you presented in your question is in a sense universal. However, the formulation $\mathbf{P_1}$ is the key to answering your question.
Suppose that $\Phi$ is an immersion of $M$ in some Euclidean space $E$ of dimension $N$ with inner product $g$. Consider the varieties
$$ N_\Phi:=\bigl\{ (x,p)\in E\times M;\;\; |v|=1,\;\;v\perp D_p\Phi(T_pM)\;\bigr\}, $$
$$ T_\Phi:=\bigl\{ (w,p)\in E\times M;\;\; |w|=1,\;\;w\in D_p\Phi(T_pM)\;\bigr\},$$
Note that $\dim N_\Phi= N-1$. Denote by $S(E)$ the unit sphere in $E$. We have two maps
$$F_\Phi: N_\Phi\to S(E),\;\;(v,p)\mapsto v, $$
$$G_\Phi: T_\Phi\to S(E),\;\; (w,p)\mapsto w. $$
Problem $\mathbf{P_1}$ is equivalent to finding two immersions $\Phi,\Psi: M\to E$ such that the ranges of $F_\Phi$ and $G_\Psi$ are disjoint. We will show that this impossible on acount of the fact that $F_\Phi$ is onto. To see this choose $v\in S(E)$. The linear map $L_v: E\to \bR$, $L_v(x)= g(v,x)$ ($g$=the inner product on $E$) induces the function $L_v^\Phi=L_v\circ \Phi$ on $S^1$. If $p\in M$ is a critical point of $L_v^\Phi$, then $(v,p)\in N_\Phi$ which shows that $v$ is in the range of $F_\Phi$.
Remark. Suppose that we require a weaker nondegeneracy condition, namely that the bilinear from $B_f(x,y)$ be nondegenerate when $x=y$. Then, as explained in this post, this bilinear form defines a canonical torsion free connection $\nabla^f$ on $TM$. More precisely, if we set $\newcommand{\pa}{\partial}$
$$ s_{jk}(x,y)=\pa^2_{x^jy^k}f(x,y), $$
and we denote by $(s^{k\ell}(x))$ the inverse of $(s_{k\ell}(x))$, then the Christoffel symbols of this connection are given by
$$ \Gamma^j_{ki}(x)=\sum_{\ell} s^{j\ell}(x)\frac{\pa^3}{\pa x^\ell\pa y^k\pa y^i}f(x,x).
$$
The question is then for which manifolds $M$ there exist smooth functions $f:M\times M\to \bR$ such that, for any $x\in M$ the matrix $$\bigl(\; \pa^2_{x^jy^k}f(x,x)\;\bigr)_{1\leq j,k\leq m}$$
is invertible.
Here is a universal example. Suppose that $M$ is a submanifold of $\bR^N$. Then the function $F:\bR^n\to\bR^n\to\bR$, $F(x,y)=x\cdot y$ restricts to a function $f$ on $M\times M$ satisfying the above weaker nondegeneracy condition. In fact, for $x\in M$ the bilinear form $B_f(x,x)$ on $T_xM$ is the inner product on $T_xM$ induced by the Euclidean inner product on $\bR^N$ so that $B_f(x,x)$ is the induced Riemann metric. As I explained in this paper the above connection is compatible with the metric and, since it is also torsion free, it must be the Levi-Civita connection !
