Is there any proof for the CLT not using characteristic functions, a simpler method?
Maybe Tikhomirov or Stein's methods?
Something self-contained you can explain to a university student (first year of mathematics or physics) and takes less than one page?
You can prove it with Stein's method, however it's debatable if the proof is elementary. The plus side of Stein's method is you get a slightly weaker form of Berry Esseen bounds essentially for free. Also, Stein's method is nothing short of black magic! You can find an exposition of the proof in section 6 of this link. You'll find other proofs of the CLT in the link as well.
Here's a brief outline:
1) Prove, using simple integration by parts and the normal distribution density, that $Ef'(A)-Xf(A)=0$ for all continuously differentiable iff $A$ is $N(0,1)$ distributed. It's easier to show $A$ normal implies the result and a bit harder to show the converse, but perhaps it can be taken on faith.
2) More generally, if $Ef(X_n)-X_nf(X_n)\rightarrow 0$ for every continuously differentiable $f$ with $f,f'$ bounded, then $X_n$ converges to $N(0,1)$ in distribution. The proof here is again by integration by parts, with some tricks. Specifically, we need to know that convergence in distribution is equivalent to $Eg(X_n)\rightarrow E g(A)$ for all bounded continuous functions $g$. Fixing $g$, this is used to reformulate:
$$Eg(X_n)-Eg(A)=Ef'(X_n)-X_nf(X_n),$$
where one solves for $f$ using basic ODE theory, and then shows $f$ is nice. Thus if we can find such a nice $f$, by assumption the r.h.s. goes to 0, and therefore so does the left side.
3) Finally, prove the central limit theorem for $Y_n:=\frac{X_1+\cdots+X_n}{\sqrt{n}}$ where $X_i$ are iid with mean 0 and variance 1. This again exploits the trick in step 2, where for every $g$ we find an $f$ such that:
$$Eg(X_n)-Eg(A)=Ef'(X_n)-X_nf(X_n).$$
Here's how I would do it if I were in high school.
Take any probability distribution with density $f(x)$, get its mean and variance $\mu_x,\sigma_x^2$. Next, approximate it with the random variable $z$ which has the following form: $$z=\mu_x-\sigma_x+2\sigma_x\xi,$$ where $\xi$ is Bernoulli random variable with parameter $p=1/2$. You can see that $\mu_z=\mu_x$ and $\sigma_z^2=\sigma_x^2$.
Now we can look at the sum $$S_n=\sum_{i=1}^n z_i$$ $$=n(\mu_x-\sigma_x)+2\sigma_x\sum_{i=1}^n\xi_i$$
You can recognize the Binomial distribution here: $\eta=\sum_{i=1}^n\xi_i$, where $\eta\sim B(n,1/2)$. You don't need characteristic function to see that it converges to normal distribution's shape.
So, in some regard you could say that the Bernoulli is the least precise approximation for any distribution, and even it converges to normal.
For instance, you can show that the moments match normal. Let's define look at the variable: $y=(S_n/n-\mu_x)\sqrt n$
$$y=\sigma_x(-1+2\eta/n)\sqrt n$$
Let's see what's the mean and variance: $$\mu_y=\sigma_x(-1+2(n/2)/n)\sqrt n=0$$ $$Var[y]=\sigma_x^2Var[2\eta/n] n=4\sigma_x^2/nn(1/4)=\sigma_x^2$$
The skewness and excess kurtosis converge to zero with $n\to\infty$, it's easy to show by plugging the known formulae for Binomial.
