Update: I added a previous version (Version 2) which does not use dimension theory and tubular neighborhoods and it might be possible to generalize it to supermanifolds. However, it requires that $N$ can be covered by finitely many coordinate balls.
VERSION 1: FOR SMOOTH MANIFOLDS AND ANY N AND $\psi$ (uses that $N$ can be covered by $\mathrm{dim}(N)+1$ charts and that there is a tubular neighborhood of $N$)
It holds
$$
\mathcal{I} = \Bigl\{ \sum_{i=1}^k (g_i \circ \psi) h_i \ \Bigl|\ k\in\mathbb{N}, g_i\in C^\infty(M), h_i\in C^\infty(L): g_i(N) = 0 \Bigr\}.
$$
Given a locally finite collection of smooth functions $(f_i\in \mathcal{I} \mid i\in I)$, we will show that
$$ \sum_{i\in I} f_i \in \mathcal{I}. $$
Lemma 1: We can assume that $f_i = (g_i \circ \psi) h_i$ for all $i\in\mathcal{I}$.
Proof: Pick a partition of unity $(\eta_j \mid j\in \mathcal{J})$ on $L$ such that $\eta_j$ has compact support for every $j\in\mathcal{J}$.
Let $\mathcal{A}:= \mathcal{I}\times\mathcal{J}$, and define
$$ f_\alpha:= \eta_j f_i $$
for all $\alpha=(i,j)\in \mathcal{A}$.
The system $(f_{\alpha}\mid \alpha\in\mathcal{A})$ is locally finite, and it holds
$$
\sum_{\alpha\in\mathcal{A}} f_{\alpha} = \sum_{j\in \mathcal{J}} \eta_j \sum_{i\in\mathcal{I}} f_i = \sum_{i\in\mathcal{I}} f_i.
$$
For every $i\in\mathcal{I}$, there is an $m_i\in \mathbb{N}$ such that $f_i = \sum_{l=1}^{m_i} (g_{il}\circ\psi)h_{il}$.
It follows that for every $\alpha = (i,j)\in\mathcal{A}$, it holds
$$
f_\alpha = \eta_j f_i = \sum_{l=1}^{m_i} (g_{il}\circ\psi)(\eta_j h_{il}) = \sum_{i=1}^{m_i} (g_{il}\circ\psi)h_{\alpha l} \in \mathcal{I},
$$
where we defined
$$ h_{\alpha l} := \begin{cases} 0 & \text{if }f_\alpha=0, \\
\eta_j h_{il} & \text{otherwise.}
\end{cases} $$
Let $x\in L$.
There is an open neighborhood $U$ of $x$ and a finite subset $\mathcal{J}_0\subset \mathcal{J}$ such that $\mathrm{supp}(\eta_j) \cap U = 0$ for all $j\in \mathcal{J}\backslash\mathcal{J}_0$.
Because $\sum_{j\in \mathcal{J}}\eta_j = 1$, it holds $U\subset \bigcup_{j\in \mathcal{J}_0} \{ x\in M \mid \eta_j(x)\neq 0\}$.
Because $\bigcup_{j\in \mathcal{J}_0} \mathrm{supp}(\eta_j)$ is compact, there exists a finite subset $\mathcal{I}_0\subset \mathcal{I}$ such that
$$
\mathrm{supp}(f_i)\cap \bigcup_{j\in \mathcal{J}_0} \mathrm{supp}(\eta_j) = \emptyset
$$ for all $i\in \mathcal{I}\backslash\mathcal{I}_0$.
Suppose that $\mathrm{supp}(h_{\alpha l})\cap U \neq 0$ for some $\alpha = (i,j)$ and $l\in \{1,\dotsc,m_i\}$.
Because $\mathrm{supp}(h_{\alpha l})\subset\mathrm{supp}(\eta_j)$, it must hold $j\in \mathcal{J}_0$.
By definition of $h_{\alpha l}$, it holds either $h_{\alpha l} = 0$ if $f_\alpha = 0$, which is equivalent to $\{x\in M \mid f_i(x)\neq 0\}\cap\{x\in M\mid \eta_j(x)\neq 0\}=\emptyset$ which is equivalent to $\mathrm{supp}(f_i)\cap\{x\in M\mid \eta_j(x)\neq 0\} = \emptyset$, or $h_{\alpha_l} = \eta_j h_{il}$.
The second option can possibly occur only for $i\in\mathcal{I}_0$.
This shows that the collection
$$ ((g_{il}\circ\psi)h_{\alpha l} \mid \alpha=(i,j)\in\mathcal{A}, l\in\{1,\dotsc,m_i\}) $$
is locally finite.
Its sum equals $\sum_{\alpha\in\mathcal{A}} f_\alpha$ and hence $\sum_{i\in \mathcal{I}}f_i$ by construction. QED
By Lemma 1, we can assume that $f_i = (g_i\circ\psi)h_i$ for $g_i\in C^\infty(M)$ with $g_i(N)=0$ and $h_i\in C^\infty(L)$ without lost of generality.
Denote $k:=\dim(N)$ and $n:=\dim(M)$.
Pick a tubular neighborhood $\mathcal{N}(N)$ of $N$ in $M$.
The $k$-dimensional manifold $N$ can be always covered by $k+1$ (not necessarily connected) charts $U_1$, $\dotsc$, $U_{k+1}$.
Every chart $U_j$ on $N$ induces a submanifold chart $V_j = \mathcal{N}(U_j)$ on $M$. Let $V_0\subset M$ be an open subset disjoint from $N$ such that $M = \cup_{j=0}^{k+1} V_j$.
Let $\lambda_0$, $\dotsc$, $\lambda_{k+1}$ be a subordinate partition of unity. Let $\mu$ be a bump function which equals $1$ on $\mathrm{supp}(\lambda_0)$ and vanishes on $N$.
Let $(x_j,y_j)\in \mathbb{R}^n$ be coordinates on $\mathcal{N}(U_j)$ such that $x_j = (x_j^1,\dotsc,x_j^k)$ gives coordinates on the base and $y_j = (y_j^1,\dotsc,y_j^{n-k})$ on fibers.
An important feature of $\mathcal{N}(U_j)$ is that it contains the vertical line $\gamma(t) = (x_j,0) + t((x_j,y_j)-(x_j,0))$ connecting $(x_j,0)$ and $(x_j,y_j)$.
The Fundamental theorem of calculus in the form
$$ f(\gamma(1))-f(\gamma(0)) = \int_{0}^1 (\nabla f)(\gamma(t))\cdot\gamma'(t) \mathrm{d}t $$
then asserts that the following holds for all $i\in I$ and $j\in \{1,\dotsc,k+1\}$ on the entire $\mathcal{N}(U_j)$:
$$ (\lambda_j g_i)(x_j,y_j) - \underbrace{(\lambda_j g_i)(x_j,0)}_{=0} = \sum_{a=1}^{n-k} y^a_j \underbrace{\int_{0}^1 \frac{\partial(\lambda_j g_i)}{\partial y^a_j}(x_j,ty_j) \mathrm{d}t}_{\displaystyle=:u_{i a}^j}. $$
Let $\tilde{y}^a_j$ and $\tilde{u}_{ia}^j$ be the smooth functions on $M$ obtained from $y^a_j$ and $u_{ia}^j$, respectively, by multiplication with a bump function which is $1$ on $\mathrm{supp} \lambda_j$ and $0$ on a neighborhood of the closure of the complement of $\mathcal{N}(U_j)$.
For all $i\in I$ and $j\in \{1,\dotsc,k+1\}$, we have the following relations on $M$:
$$ \lambda_0 g_i = \mu \lambda_0 g_i\quad\text{and}\quad\lambda_j g_i = \sum_{a=1}^{n-k} \tilde{y}^a_j \tilde{u}_{ia}^j. $$
Using this, we compute
\begin{align*}
\sum_{i\in I} (g_i \circ \psi) h_i
&= \sum_{i\in I} \sum_{j=0}^{k+1} (\lambda_j g_i \circ \psi) h_i \\
& = \sum_{i\in I} (\lambda_0 g_i \circ \psi) h_i + \sum_{i\in I} \sum_{j=1}^{k+1} \sum_{a=1}^{n-k} (\tilde{y}^a_j \circ \psi)(\tilde{u}_{ia}^j \circ \psi)h_i \\
& = (\mu\circ\psi)\sum_{i\in I}(\lambda_0 g_i \circ \psi) h_i+ \sum_{j=1}^{k+1} \sum_{a=1}^{n-k} (\tilde{y}^a_j\circ \psi) \sum_{i\in I} (\tilde{u}_{ia}^j \circ \psi)h_i\\
& = (G_0 \circ \psi) H_0 + \sum_{j=1}^{k+1} \sum_{a=1}^{n-k} (G_{ja}\circ\psi)H_{ja},
\end{align*}
where we denoted
$$ G_0:= \mu,\quad G_{ja}:=\tilde{y}^a_j,\quad H_0:=\sum_{i\in I}(\lambda_0 g_i \circ \psi) h_i,\quad H_{ja}:=\sum_{i\in I} (\tilde{u}_{ia}^j\circ \psi)h_i. $$
It holds $G_0$, $G_{ja}\in C^\infty(M)$, $G_0(N)=G_{ja}(N) = 0$, $H_0$, $H_{ja}\in C^\infty(L)$, and it follows that $\sum_{i\in I} (g_i \circ \psi) h_i \in \mathcal{I}$.
VERSION 2: PROOF WHEN $N$ CAN BE COVERED BY FINITELY MANY COMPATIBLE COORDINATE BALLS (not using dimension theory and tubular neighborhood)
Write $\mathbb{R}^n = \mathbb{R}^k\times \mathbb{R}^{n-k}$ with coordinates $(x,y)$. Let $f: \mathbb{R}^n\rightarrow \mathbb{R}$ be a smooth function vanishing at $\{(x,y) \mid x = 0\}$. Then the fundamental theorem of calculus asserts that the following holds for all $(x,y)\in \mathbb{R}^n$: $$ f(x,y) = \sum_{j=1}^{k} x^j \int_{0}^1 \frac{\partial f}{\partial x^j}(tx,y) dt. $$ Let $U_\alpha$ $(\alpha\in\mathcal{A})$ be a cover of $N$ by coordinate balls, and let $\lambda_\alpha$ $(\alpha\in\mathcal{A})$ be a subordinate partition of unity. Suppose that we are given $\sum_{i\in I} (g_i \circ \psi) h_i$ as above and we want to show that it lies in $\mathcal{I}$. We can even assume that the support lies in an arbitrary small neighborhood of $\psi^{-1}(N)$. Using the analytical fact above, there are smooth functions $x_{\alpha}^j$ vanishing on $N$ and smooth functions $u^{\alpha}_{ij}$ for all $i\in I$, $\alpha\in\mathcal{A}$ and $j\in\{1,\dotsc,k:=\mathrm{codim} N\}$ such that $$ \lambda_\alpha g_i = \sum_{j=1}^k x_\alpha^j u_{ij}^\alpha. $$ We compute \begin{align*} \sum_{i\in I} (g_i \circ \psi) h_i = \sum_{j=1}^k \sum_{\alpha\in\mathcal{A}} (\lambda_\alpha x^j_\alpha\circ\psi) \sum_{i\in I} (u_{ij}^\alpha \circ \psi) h_i. \end{align*} If $\mathcal{A}$ is finite, then we are done.