$\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$
