[proofplan]
The argument is the relative version of Moser's method. First we use the ambient agreement of $\omega_0$ and $\omega_1$ along $N$ to write $\omega_1-\omega_0=d\sigma$ with a primitive $\sigma$ that vanishes as a covector on every tangent space $T_xM$ for $x\in N$. Then we solve the time-dependent equation $\iota_{X_t}\omega_t=-\sigma$; the vanishing of $\sigma$ along $N$ forces $X_t$ to vanish along $N$, so its flow fixes $N$ pointwise. Finally, after shrinking neighbourhoods so the flow is defined up to time $1$, the usual Moser derivative computation shows that the time-one map pulls $\omega_1$ back to $\omega_0$.
[/proofplan]
[step:Produce a relative primitive that vanishes as an ambient covector along $N$]
Let $W\subset M$ be an open neighbourhood of $N$ on which $\omega_0$, $\omega_1$, and every $\omega_t$ are defined and symplectic. Define the smooth $2$-form
\begin{align*}
\alpha:=\omega_1-\omega_0\in\Omega^2(W).
\end{align*}
Since $\omega_0$ and $\omega_1$ are symplectic, they are closed, and therefore
\begin{align*}
d\alpha=d\omega_1-d\omega_0=0.
\end{align*}
The hypothesis gives $\alpha_x=0$ as an alternating form on $T_xM$ for every $x\in N$.
Since $N\subset M$ is an embedded submanifold, choose a tubular neighbourhood of $N$ in $M$ and shrink $W$ inside it. The relative homotopy operator used in the [Relative Poincare Lemma][citetheorem:10050] is computed along the fibres of the tubular projection onto $N$. Applied to the closed form $\alpha$ with the stronger ambient vanishing condition $\alpha_x=0$ on $T_xM$, it gives a smooth $1$-form
\begin{align*}
\sigma\in\Omega^1(W)
\end{align*}
such that
\begin{align*}
d\sigma=\alpha.
\end{align*}
Moreover, at a point $x\in N$, the homotopy integral contains only evaluations of $\alpha_x$ with arguments in $T_xM$ and at least one fibre-radial vector; since $\alpha_x=0$ as a [bilinear form](/page/Bilinear%20Form) on all of $T_xM$, the resulting covector satisfies
\begin{align*}
\sigma_x=0
\end{align*}
as a covector on $T_xM$ for every $x\in N$. We henceforth replace $W$ by this smaller neighbourhood.
[/step]
[step:Solve the Moser equation and prove the vector field vanishes on $N$]
For each $t\in[0,1]$ and each $x\in W$, nondegeneracy of $(\omega_t)_x$ gives an isomorphism $\Lambda_{t,x}:T_xM\to T_x^*M$ defined by $\Lambda_{t,x}(v)=\iota_v(\omega_t)_x$.
Define the smooth time-dependent vector field $X:[0,1]\times W\to TW$, written $X(t,x)=X_t(x)$, by the equation
\begin{align*}
\iota_{X_t}\omega_t=-\sigma.
\end{align*}
Smoothness follows because the inverse of the smooth bundle isomorphism induced by $\omega_t$ depends smoothly on $(t,x)$.
If $x\in N$, then $\sigma_x=0$, so
\begin{align*}
\iota_{X_t(x)}(\omega_t)_x=0.
\end{align*}
Since $(\omega_t)_x$ is nondegenerate, this implies
\begin{align*}
X_t(x)=0
\end{align*}
for every $t\in[0,1]$ and every $x\in N$.
[guided]
The equation
\begin{align*}
\iota_{X_t}\omega_t=-\sigma
\end{align*}
is the central Moser equation. For fixed $t\in[0,1]$ and $x\in W$, the form $(\omega_t)_x$ is symplectic, hence nondegenerate. Nondegeneracy means that the [linear map](/page/Linear%20Map) $\Lambda_{t,x}:T_xM\to T_x^*M$ defined by $\Lambda_{t,x}(v)=\iota_v(\omega_t)_x$ is an isomorphism. Therefore there is a unique vector $X_t(x)\in T_xM$ satisfying
\begin{align*}
\iota_{X_t(x)}(\omega_t)_x=-\sigma_x.
\end{align*}
As $(t,x)$ varies, the coefficients of $\omega_t$ and $\sigma$ vary smoothly, and inversion of an invertible matrix is smooth; hence the resulting time-dependent vector field $X$ is smooth.
The point of requiring the primitive to vanish as an ambient covector along $N$ now appears. If $x\in N$, then $\sigma_x=0$ on all of $T_xM$, not merely on $T_xN$. Thus the Moser equation becomes
\begin{align*}
\iota_{X_t(x)}(\omega_t)_x=0.
\end{align*}
Because $(\omega_t)_x$ is nondegenerate on the full tangent space $T_xM$, the only vector whose contraction with $(\omega_t)_x$ is zero is the zero vector. Hence
\begin{align*}
X_t(x)=0
\end{align*}
for every $x\in N$ and every $t\in[0,1]$.
[/guided]
[/step]
[step:Shrink the domain so the flow exists through time $1$ and fixes $N$]
For each $x\in N$, the constant curve $t\mapsto x$ solves the time-dependent [initial value problem](/page/Initial%20Value%20Problem) on the compact interval $[0,1]$ because $X_t(x)=0$ for all $t\in[0,1]$. By the local existence theorem with continuous dependence for smooth time-dependent vector fields, applied along this solution on the compact time interval $[0,1]$, there is an open neighbourhood $V_x\subset W$ of $x$ such that every solution starting in $V_x$ is defined for all $t\in[0,1]$. Define
\begin{align*}
U_0:=\bigcup_{x\in N}V_x.
\end{align*}
Then $U_0\subset W$ is an open neighbourhood of $N$, and the flow $\varphi:[0,1]\times U_0\to W$, written $\varphi(t,y)=\varphi_t(y)$, is defined for every $t\in[0,1]$. Moreover, for every $x\in N$, uniqueness of solutions applied to the constant solution gives
\begin{align*}
\varphi_t(x)=x
\end{align*}
for every $x\in N$ and every $t\in[0,1]$.
Set
\begin{align*}
U_1:=\varphi_1(U_0).
\end{align*}
After shrinking $U_0$ if necessary, $\varphi_1:U_0\to U_1$ is a diffeomorphism onto the open neighbourhood $U_1$ of $N$.
[/step]
[step:Compute the derivative of the pulled-back forms]
For $t\in[0,1]$, define the pulled-back $2$-form
\begin{align*}
\beta_t:=\varphi_t^*\omega_t\in\Omega^2(U_0).
\end{align*}
The time-dependent pullback differentiation formula gives
\begin{align*}
\frac{d}{dt}\beta_t=\varphi_t^*(\dot{\omega}_t+\mathcal L_{X_t}\omega_t).
\end{align*}
Here
\begin{align*}
\dot{\omega}_t=\omega_1-\omega_0=\alpha=d\sigma.
\end{align*}
By Cartan's formula for the Lie derivative and by closedness of $\omega_t$,
\begin{align*}
\mathcal L_{X_t}\omega_t=d(\iota_{X_t}\omega_t)+\iota_{X_t}d\omega_t=d(-\sigma)+0=-d\sigma.
\end{align*}
Therefore
\begin{align*}
\dot{\omega}_t+\mathcal L_{X_t}\omega_t=d\sigma-d\sigma=0.
\end{align*}
Thus
\begin{align*}
\frac{d}{dt}\beta_t=0
\end{align*}
for every $t\in[0,1]$, so $\beta_t$ is independent of $t$.
[/step]
[step:Take the time-one map to obtain the desired diffeomorphism]
Since $\beta_t=\varphi_t^*\omega_t$ is independent of $t$, evaluating at $t=0$ and $t=1$ gives
\begin{align*}
\varphi_1^*\omega_1=\varphi_0^*\omega_0.
\end{align*}
The flow satisfies $\varphi_0=\operatorname{id}_{U_0}$, so
\begin{align*}
\varphi_1^*\omega_1=\omega_0.
\end{align*}
Define
\begin{align*}
\psi:=\varphi_1:U_0\to U_1.
\end{align*}
From the previous step, $\psi$ is a diffeomorphism. Since the flow fixes $N$ pointwise,
\begin{align*}
\psi|_N=\operatorname{id}_N.
\end{align*}
Together with
\begin{align*}
\psi^*\omega_1=\omega_0,
\end{align*}
this is precisely the required conclusion.
[/step]