$\newcommand{\R}{\mathbb R}\newcommand{\la}{\lambda}\newcommand{\Ga}{\Gamma}$We have \begin{equation*} f_{j+1}(x)=(f_j*f_j)(x)\,1(\|x\|<1)/A_j \end{equation*} for $j=0,1,\dots$ and $x\in\R^n$, where $A_j:=\int_{\|x\|<1}dx\,(f_j*f_j)(x)$.
For each $j$ and each $x$, the value of $(f_j*f_j)(x)$ is an $n$-fold integral.
If $f_0$ is spherically symmetric, then so will be $f_j$ for each $j$. Moreover, then, for $n\ge2$, the $n$-fold integral can be reduced to a double one.
Indeed, in what follows, do assume that $f_0$ is spherically symmetric. Assume also that $n\ge2$. Then for each $j=0,1,\dots$ there is a measurable function $g\colon[0,\infty)\to[0,\infty)$ such that \begin{equation*} f_j(x)=g_j(\|x\|) \end{equation*} for all $x\in\R^n$. So, \begin{equation*} (f_j*f_j)(x)=(Gg)(\|x\|), \end{equation*} where \begin{equation*} (Gg)(u):=\int_0^1\int_0^1\la_n(E_{u,a,b,da,db})g(a)g(b), \end{equation*} where $\la_n$ is the Lebesgue measure over $\R^n$, \begin{equation*} E_{u,a,b,da,db}:=\{y\in\R^n\colon a<\|y\|<a+da,\,b<\|ue_1-y\|<b+db\}, \end{equation*} and $e_1:=(1,0,\dots,0)\in\R^n$.
Shaded in the picture below is the small approximate parallelogram that is the intersection of the set $E_{u,a,b,da,db}$ (with $(u,a,b,da,db)=(6,5,4,0.4,0.4)$) with the "upper" half-plane of the 2D plane through the points $0$, $ue_1$, and $y$, where $y$ is such that $\|y\|=a$ and $\|ue_1-y\|=b$; the "upper" half-plane is bounded by the line through $0,ue_1$ and contains the point $y$. The value of $r$ is the distance from the point $y$ to the line through $0,ue_1$.
For infinitesimal $da$ and $db$, the area of the parallelogram is $dA=da\,db\,\sin t$, where $t$ is the angle between the vectors $y$ and $ue_1-y$. Also, expressing the area of the triangle $0(ue_1)y$ in two different ways, we see that $r=\frac{ab}u\,\sin t$. Also, by the cosine theorem, $\cos t=\frac{a^2+b^2-u^2}{2ab}$ and hence $\sin t=\frac1{2u}\,\sqrt{4a^2b^2-(a^2+b^2-u^2)^2}$ for $u>0$. We also have the triangle inequalities $|a-b|\le u\le a+b$. So, \begin{equation*} \begin{aligned} &\la_n(E_{u,a,b,da,db})\\ &=dA\,c_{n-2}r^{n-2}\,1(|a-b|\le u\le a+b) \\ &=\frac{c_{n-2}}{2^{n-1}}\frac{da\,db}{abu^{n-2}}\,(4a^2b^2-(a^2+b^2-u^2)^2)^{(n-1)/2} \,1(|a-b|\le u\le a+b), \end{aligned} \end{equation*} where $c_k=\dfrac{2\pi^{(k+1)/2}}{\Ga((k+1)/2)}$ is the "length/surface area" of the $k$-dimensional unit sphere $\mathbb S_k$ in $\R^{k+1}$, so that $c_0=2$, $c_1=2\pi$, etc. So, for a certain positive real constant $b_n$ depending only on $n$ and all $u\in(0,1]$, \begin{equation*} (Gg)(u)=\frac{b_n}{u^{n-2}} \int_0^1\frac{da\,g(a)}a\, \\ \times\int_{\max(a,u-a)}^{\min(1,u+a)} \frac{db\,g(b)}b\,(4a^2b^2-(a^2+b^2-u^2)^2)^{(n-1)/2}. \end{equation*}
Finally, for $j=0,1,\dots$ and $u\in[0,1]$, \begin{equation*} g_{j+1}(u)=(Hg_j)(u), \end{equation*} where \begin{equation*} (Hg)(u):=(Gg)(u)\Big/A(g) \end{equation*} and \begin{equation*} A(g):=\int_{\|x\|<1}dx\,(Gg)(\|x\|)=c_{n-1}\int_0^1 du\,u^{n-1}(Gg)(u). \end{equation*}
The case $n=1$ is exceptional and easier.
Below are the graphs $\{(u,g_j(u))\colon0<u\le1\}$ for constant $g_0$, for $j=1$ (red), $j=2$ (green), $j=3$ (blue), and for $n\in\{1,2,3,5,10,20,100\}$ (the amount of calculations is about the same for all these $n$):
It appears that the $g_j$'s converge very fast (strangely enough, the fastest convergence seems to be for $n=5$, as compared with the other $n\in\{1,2,3,5,10,20,100\}$).