Implementable numerical scheme for the equation $a=\text{Erf}\big(z/\sqrt{2N_{a}}\big)$ Let $z>0$ be fixed and $A$ be the set of non-increasing functions from $\mathbb R_+$ to $[0,1]$ with norm $\|\cdot\|:=\|\cdot\|_\infty$. Define by $F$ the operator on $A$ by
\begin{equation*}
    F(a)(t):=\text{Erf}\left(\frac{z}{\sqrt{2N_a(t)}} \right),\quad \forall t\ge 0,  
\end{equation*}
where $\text{Erf}$ is the Gauss error function and $N_a:\mathbb R_+\to\mathbb R_+$ is defined by
$$N_a(t):=\int_0^t\frac{ds}{(1+a(s))^2}.$$
Iosif has shown in Numerical solution to some functional equation that $F$ maps $A$ to $A$ and is a contraction map. I'm looking for a numerical approximation of the fixed point $a^*$ of $F$.  Let $a_0\equiv 0$ and $a_{n}:=F(a_{n-1})$ for all $n\ge 1$. Then it follows that
$$\|a_n-a^*\|\le Cr^n,\quad \forall n\ge 1,$$
where $C>0, r\in (0,1)$ are some constants. How could we implement (code in computer) the iteration $a_{n}:=F(a_{n-1})$, and obtain a numerical approximation of $a_N$ for any $N$?
PS : Here the numerical approximation of $a_N$ means that: For any fixed $T>0$, any $\epsilon>0$ and any subset $\{0=t_0<t_1<\cdots<t_M=T\}$, the computer yields a output $\{y_0, y_1,\ldots, y_M\}$ s.t.
$$\max_k |y_k-a_N(t_k)|\le \epsilon.$$
Of course, any other numerical method (different from the above iteration) is highly appreciated.
 A: $\newcommand\erf{\operatorname{erf}}\newcommand\R{\mathbb R}\newcommand{\de}{\delta}\newcommand\ep\epsilon
$In the previous answer, it was shown that the operator $F$ on $A$ is $r$-Lipschitz for a certain universal constant $r\in(0,1)$ with respect to the norm $\|\cdot\|_\infty$. The same proof holds for the corresponding operator $F_T\colon A_T\to A_T$, where $A_T$ is the set of all non-increasing functions from $(0,T]$ to $[0,1]$ with norm given by $\|a\|:=\sup_{t\in(0,T]}|a(t)|$ for $a\in A_T$ and
\begin{equation*}
    F_T(a)(t):=\text{Erf}\left(\frac{z}{\sqrt{2N_a(t)}} \right)   
\end{equation*}
for $a\in A_T$ and $t\in(0,T]$; this follows because for any $t\in(0,T]$ the value of $N_a(t)$ depends only on the values of $a$ on $(0,T]$.
The additional ingredient needed to answer the current question is the observation that $F_T(a)(t)$ is $L$-Lipschitz in $t$ with the same $L$ for all $a\in A_T$. Indeed, for any $a\in A_T$ and $t\in(0,T]$, with $u:=z/\sqrt{N_a(t)}$,
\begin{equation*}
    F_T(a)'(t):=\frac d{dt}\,F_T(a)(t)
    =-\frac2{\sqrt\pi}\,e^{-u^2/2}u^3\,\frac1{2\sqrt2\,z^2}\frac1{(1+a(t))^2},
\end{equation*}
so that
\begin{equation*}
    |F_T(a)'(t)|\le L:=c/{z^2},
\end{equation*}
where $c:=(3/e)^{3/2}/\sqrt{2\pi}$.
As requested, take now any real $\ep>0$, any natural $K$, and any $t_0,\dots,t_K$ such that  $0=t_0<t_1<\cdots<t_K=T$. For each $n=0,1,\dots$, we want a computer to yield $y_{n,0},\ldots,y_{n,K}$ s.t.
\begin{equation*}
    M_n:=\max_{k\in[K]}|a_n(t_k)-y_{n,k}|\le\ep, \tag{$*$}\label{1}
\end{equation*}
where, as usual, $[K]:=\{1,\dots,K\}$.
Take any real $\de>0$. By refining the "partition" $0=t_0<t_1<\cdots<t_K=T$, assume without loss of generality that
\begin{equation*}
    \max_{k\in[K]}|t_k-t_{k-1}|\le\de.
\end{equation*}
Let (say) $a_0:=1$ and $a_n:=F_T(a_{n-1})$ for $n\ge1$.
For $k\in[K]$, let $y_{0,k}:=0$ and for $n\ge1$ let (a computer compute)
\begin{equation*}
    y_{n,k}:=F_T(Y_{n-1})(t_k),
\end{equation*}
where
\begin{equation*}
    Y_n(t):=\sum_{k\in[K]}y_{n,k}\,1(t_{k-1}<t\le t_k). 
\end{equation*}
Then, recalling that $F_T$  is $r$-Lipschitz, for $n\ge1$ we have
\begin{equation*}
\begin{aligned}
M_n&=\max_{k\in[K]}|F_T(a_{n-1})(t_k)-F_T(Y_{n-1})(t_k)| \\ 
&\le\|F_T(a_{n-1})-F_T(Y_{n-1})\| \\ 
&\le r\|a_{n-1}-Y_{n-1}\| \\ 
&=r\sup_{t\in(0,T]}|a_{n-1}(t)-Y_{n-1}(t)| \\ 
&=r\max_{k\in[K]}\sup_{t\in(t_{k-1},t_k]}|a_{n-1}(t)-y_{n-1,k}|;  
\end{aligned}
\end{equation*}
also, $|a_{n-1}(t)-y_{n-1,k}|\le|a_{n-1}(t_k)-y_{n-1,k}|+L\de$, because ($F_T(a)(t)$ is $L$-Lipschitz in $t$ and hence) $a_{n-1}$ is $L$-Lipschitz.
It follows that for all $n\ge1$
\begin{equation*}
    M_n\le rM_{n-1}+rL\de
\end{equation*}
and hence (by induction on $n$, with $M_0=0$)
\begin{equation*}
    M_n\le C\de,
\end{equation*}
where $C:=rL/(1-r)$. Finally, choosing $\de=\ep/C$, we get \eqref{1}. $\quad\Box$
