I am trying to prove what I stated in the title. I am not sure I can find a counterexample either but I can not figure out how prove this directly.
So far I've tried finding an open neighbourhood in $TM$ such that $dF|_U:U \rightarrow TN$ is a smooth embedding but I do not think it is the right way. Any indication?
Thanks!
Let me write the map $dF$ in local coordinates. I will use $x_1, \dots x_n$ as the manifold coordinates around $p$ and $u_1, \dots , u_n$ as the tangent vector coordinates. I will use $y_1, \dots y_n$ as the manifold coordinates around $F(p)$ and $v_1, \dots , v_n$ as the tangent vector coordinates.
Suppose that the map $F$ is given by $$ (x_1, \dots, x_n) \mapsto (y_1, \dots , y_n) = (F_1(x_1, \dots, x_n), \dots, F_n(x_1, \dots, x_n)).$$
Then the associated map on the total space of the tangent bundle is $$(x_1, \dots, x_n,u_1, \dots, u_n) \mapsto (y_1, \dots , y_n,v_1, \dots, v_n) = \\ \left(F_1(x_1, \dots, x_n), \dots, F_n(x_1, \dots , x_n),\sum_{i}\frac{\partial F_1(x_1, \dots, x_n)}{\partial x_i} u_i \ , \ \dots \ , \sum_{i}\frac{\partial F_n(x_1, \dots, x_n)}{\partial x_i} u_i \right)$$ [Sorry, that was meant to fit onto one line.]
You can view this as a map from an open subset of $\mathbb R^n \times \mathbb R^n$ to $\mathbb R^n \times \mathbb R^n$.
Now let us find the derivative of this map. This derivative should look like a $2n \times 2n$ matrix of smooth functions: $$ \left( \begin{array}{cc} \frac{\partial F_j}{\partial x_i} & 0 \\ \sum_k \frac{\partial^2 F_j}{\partial x_i \partial x_k} u_k & \frac{\partial F_j}{\partial x_i} \end{array}\right)$$
Now suppose $F$ is a submersion at $p$. Since $M$ and $N$ have the same dimension, this implies that $dF_p$ is bijective. If the point $p$ has coordinates $(x_1, \dots , x_n) = (0, \dots, 0)$, then this means that the $n \times n$ matrix $$ \left( \frac{\partial F_j(0,\dots,0)}{\partial x_i} \right)_{ij}$$ has non-vanishing determinant.
However, the derivative of the above map from the open subset of $\mathbb R^n \times \mathbb R^n$ to $\mathbb R^n \times \mathbb R^n$, when viewed as a $2n \times 2n$ matrix of smooth functions, has determinant equal to $$ \left( \det \left( \frac{\partial F_j}{\partial x_i} \right)_{ij} \right)^2$$ When evaluated at $(x_1, \dots, x_n, u_1, \dots, u_n) = (0,\dots, 0, u_1, \dots, u_n)$ for any $(u_1, \dots, u_n)$, this determinant is certainly non-zero.
Applying the inverse function theorem, we conclude that our map $dF$ between the total spaces of the two tangent bundles is a local diffeomorphism around $(0,\dots, 0, u_1, \dots, u_n)$.