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* !