$\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_j)(\|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$. 

[![enter image description here][1]][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{2ab}\,\sqrt{4a^2b^2-(a^2+b^2-u^2)^2}$. 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$): 

[![enter image description here][2]][2]


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\}$). 


  [1]: https://i.sstatic.net/XtH32.png
  [2]: https://i.sstatic.net/259T6.png