Solution 1:

Pick coordinate charts around $p \in M$ and $f(p) \in N$; assume that $M$ and $N$ are of dimension $m$, $n$ respectively, and $($by continuity$)$ that the points of $M$ corresponding to our chart around $p$ maps into those of our chart around $f(p)$, i.e. we may reduce to the case in which $M$, $N$ are subsets of Euclidean space, and even assume that $M$, $N$ are open rectangles $($products of intervals$)$. Then, if $\pi_i$ is the projection onto the $i$th coordinate, we may pick $1 \le i_1 < i_2 < \dots < i_m \le n$ such that $\pi_{i_1} \circ f$, $\pi_{i_m} \circ f$ have linearly independent gradients $($since $df$$)$ has rank $m$. Let $\pi$ be the projection onto the subspace of $\mathbb{R}^n$ on which all but these $m$ coordinates vanish. Now, by the Inverse Function Theorem, there is a neighborhood $U$ of $p$ such that the $m$-tuple of functions $\pi_{i_j} \circ f$, $1 \le j \le m$ restricted to $U$ gives a diffeomorphism onto an open set in $\mathbb{R}^m$. Let $g$ be the inverse diffeomorphism, which is defined on an open set near $\pi(f(p))$. We then have that $g \circ \pi$ is smooth and $(g \circ \pi) \circ f$ is the identity on $U$, so $f|_U$ has a continuous inverse defined on $f(U)$. Thus $f|_U$ is a smooth homeomorphism onto its image with injective derivative, hence an embedding.

$($The link Immersions are Locally Embeddings might also be helpful.$)$