$\newcommand{\de}{\delta}
\newcommand{\De}{\Delta}
\newcommand{\ep}{\varepsilon}
\newcommand{\ga}{\gamma}
\newcommand{\Ga}{\Gamma}
\newcommand{\la}{\lambda}
\newcommand{\Si}{\Sigma}
\newcommand{\thh}{\theta}
\newcommand{\R}{\mathbb{R}}
\newcommand{\E}{\operatorname{\mathsf E}}
\newcommand{\PP}{\operatorname{\mathsf P}}
\newcommand{\ii}[1]{\operatorname{\mathsf I}\{#1\}}$

It is straightforward, but a bit tedious, to construct a function $f\in C^\infty(\R)$ satisfying conditions 1--4 and such that $f(x)=e^{-1/x}$ for $x$ in a right open neighborhood $N$ of $0$ (see details at the end of this answer). Then
\begin{equation}
F(x,y)=e^{1/z-1/x},
\end{equation}
where $z:=x+ye^{-1/x}$; everywhere here,
$y>-1$ and $x$ is a small enough positive real number (depending on $y$) such that $x$ and $z$ are both in $N$. Note that $z'_x=1+ye^{-1/x}/x^2$ and $z'_y=e^{-1/x}$, whence we have the crucial observation:
\begin{multline*}
\Big(\frac1z-\frac1x\Big)'_x
=\frac1{x^2}-\frac1{z^2}\Big(1+e^{-1/x}\frac y{x^2}\Big)
=\frac{(z-x)(z+x)}{x^2 z^2}-e^{-1/x}\frac y{x^2 z^2} \\
=y e^{-1/x}\Big(\frac1{xz^2}+\frac1{zx^2}-\frac1{x^2z^2}\Big).
\end{multline*}

So, we see that
\begin{align*}
F'_x(x,y)&=F(x,y)p_1(1/x,1/z,y,e^{-1/x})e^{-1/x},\\
F'_y(x,y)&=F(x,y)q_1(1/x,1/z,y,e^{-1/x})e^{-1/x}
\end{align*}
for some polynomials $p_1$ and $q_1$. So, by induction, all partial derivatives of $F(x,y)$ of orders $\ge1$ are of the form $F(x,y)p(1/x,1/z,y,e^{-1/x})e^{-1/x}$ for some polynomials $p$.

Also, $|z-x|\ll e^{-1/x}$, $z\sim x$, $F(x,y)\to1$, which implies that all partial derivatives of $F(x,y)$ of orders $\ge1$ go to $0$; the convergence here is for $x\downarrow0$ uniformly over all $y$ in any compact subset of $(-1,\infty)$. This implies that $F\in C^\infty(\R\times(-1,\infty))$, with $F(x,y):=1$ for $x\le0$.

**Added:** details concerning the first sentence of this answer. For all real $x$,
let
\begin{equation}
f(x) :=\int_\R g(x + \ep u(x)t)Cu(t)\,dt,
\end{equation}
where
\begin{equation}
g(x) := e^{-1/x} \ii{0 <x\le2}+[e^{-1/2} +(x-2)(1-e^{-1/2})]\,\ii{2 < x \le 3}+\ii{x > 3}
\end{equation}
\begin{equation}
u(x) := \exp\Big\{-\frac1{(x - 1)(4 - x)}\Big\}\,\ii{1 < x < 4},
\end{equation}
$C := 1/\int_\R u(x)\,dx$, $\ii{\cdot}$ is the indicator, and $\ep$ is a positive real number small
enough so that $1+4\ep u'> 0$ (whence $x+\ep u(x)t$ is increasing in $x \in\R$ for each $t$ in
the interval $(1,4)$).

Then $f$ is in $C^\infty(\R)$, satisfies conditions 1--4, and $f(x)=e^{-1/x}$ for $x\in(0,1)$. In particular, $f$ is increasing -- because $g$ is so and $x+\ep u(x)t$ is increasing in $x \in\R$ for each $t\in(1,4)$. Also, $f(x)=e^{-1/x}<x$ for $x\in(0,1)$, $g\le1$ and hence $f\le1$ on $\R$, and so, $f(x)<x$ for $x>1$. That $f$ is $C^\infty$ on the interval $(1,4)$ follows because
\begin{equation}
f(x) =\int_\R g(y)Cu\Big(\frac{y-x}{\ep u(x)}\Big)\,\frac{dy}{\ep u(x)}
\end{equation}
for $x\in(1,4)$.