I'm trying to understand how individual contributions of single particles can add up to something one would call a phase transition. For this I consider a model of a system of $N$ one-dimensional particles in a potential well like this:
$$U(x_i)=\begin{cases} 0&\text{if }|x|<d\\ U_{II}&\text{if }d<|x|<L\\ \infty&\text{if }|x|>L \end{cases}.$$
Here $x_i$ is position of $i$th particle, $U_{II}>0$ is height of potential barrier, $d$ is width of potential well and $L>d$ is the width of the total motion area. There's no interaction between the particles, only with the external potential $U$.
Now I calculate how average kinetic energy $\left\langle T\right\rangle$ depends on average total energy per particle. This should be proportional to temperature, since both inside the well and outside of it I have an ideal or almost ideal gas. Here's what I get for $\langle T\rangle$:
$$\langle T\rangle=E-U_{II}\frac{\Delta t_{II}}{\Delta t_I+\Delta t_{II}}=\begin{cases}\frac{(L-d)U_{II}}{L+d\left(1-\sqrt{1-\frac{U_{II}}{E}}\right)}+E&\text{if }E>U_{II}\\ E&\text{if }E<U_{II}\end{cases},$$
where $\Delta t_{I,II}$ are times spent in the areas of $U=0$ and $U=U_{II}$ respectively per period of motion. Here's what it looks like for $N=1,2,3,4,5,10$ (yes, the last one is $10$ not $6$), $L=30$, $d=0.1$ and $U_{II}=7$:
I calculate it as
$$\langle T\rangle=\int_0^\infty dE_1\int_0^\infty dE_2\cdots\int_0^\infty dE_N\left(\frac1N\sum_{i=1}^N\langle T_i\rangle\delta(NE-E_1-E_2-\cdots-E_N)\right)$$
using Monte-Carlo integration. Here $\langle T_i\rangle$ is kinetic energy of $i$th particle averaged in time, and $\langle T\rangle$ is $\langle T_i\rangle$ averaged over all particles and possible combinations of energies of each particle $E_i$ such that total energy of the system is $EN$. The discontinuity of the curve for one particle is due to the particle going out of the potential well ($E\gtrsim U_{II}$): near this energy the particle most of the time moves very slowly above the barrier of $U=U_{II}$, thus its average kinetic energy is small.
And here's how what I suppose to be heat capacity, namely $\frac{dE}{d\langle T\rangle}$, looks (the noise is due to derivatives of Monte-Carlo-integrated function)
This seems to resemble the peak of heat capacity of a second order transition, c.f. these plots for the Ising model simulations.
But my doubt is that as I increase number of particles, the peak appears to decrease its amplitude and widen. This shouldn't happen for a phase transition, which is normally a many-particle effect. So my question is: is what I described above even a phase transition at all? Or will the peak of heat capacity actually drop to zero once I take the limit of $N\to\infty$? Or am I completely wrong about what I think is heat capacity or temperature for this system?