Let $G$ be a Lie group and $l_g:G\rightarrow G$ be the left-multiplication map.
Let $X$ be a left invariant vector field on $G$, i.e., $X:G\rightarrow TG$ is such that $(l_g)_*X=X$ on $G$ where $(l_g)_*X$ is defined as follows: $$((l_g)_*X)_{gh}=(l_g)_{*,h}(X_h)$$ where $(l_g)_{*,h}:T_hG\rightarrow T_{gh}G$ is the differential of $l_g$ at $h$.
So, supposing $X$ is a left-invariant vector field, we have $$(l_g)_{*,e}(X_e)=X_g.$$
I am trying to prove that the Lie bracket $[X,Y]$ is also a left-invariant vector field.
$${[X,Y]_g(f)}={X_g(Y(f))-Y_g(X(f))}={(l_g)_{*,e}(X_e)(Y(f))-(l_g)_{*,e}(Y_e)(X(f)) }={X_e(Y(f)\circ l_g)-Y_e(X(f)\circ l_g)}={X_e(Y(f\circ l_g))-Y_e(X(f\circ l_g))}={[X,Y]_e(f\circ l_g)}={(l_g)_{*,e}([X,Y]_e)(f)}$$
This is true for all $f$. So, we have $$(l_g)_{*,e}([X,Y]_e)=[X,Y]_g.$$
Using same idea, we can prove that $$(l_g)_{*,h}([X,Y]_h)=[X,Y]_{gh}$$ for all $h\in G$ which is same as saying that $[X,Y]$ is a left-invariant vector field.
I am sure about all equalities except one, i.e., $Y(f\circ l_g)=Y(f)\circ l_g$.
I am sure this is correct but could not see. Any suggestions about the proof are welcome.
I have not gone through all your computations, I did just check the equality you were doubting.
First, let me recall that by definition, for all $x\in G$ and all smooth map $h\colon G\rightarrow\mathbb{R}$, one has the following: $$Y(h)(x):=T_xh(Y(x))\tag{1}.$$ Therefore, using the definition $(1)$ and the chain rule formula, one gets: $$Y(f\circ l_g)(x)=T_x(f\circ l_g)(Y(x))=T_{gx}f(T_xl_g(Y(x)))\tag{2}.$$ Besides, using twice the left-invariance of $Y$ and the chain rule formula along with $l_{gx}=l_g\circ l_x$, one has: $$Y(gx)=T_el_{gx}(Y(e))=T_e(l_g\circ l_x)(Y(e))=T_xl_g(T_el_x(Y(e))=T_xl_g(Y(x))\tag{3}.$$ Hence, using the definition $(1)$ and the equality $(3)$, one gets: $$Y(f)(l_g(x))=T_{gx}f(Y(gx))=T_{gx}f(T_xl_g(Y(x)))\tag{4}.$$ Finally, the left terms of $(2)$ and $(4)$ are equal. Whence the result.
Here is another way to see that if $X$ and $Y$ are left-invariant vector fields of $G$, then $[X,Y]$ is left invariant.
This identity is easily checked on the level of derivations of $C^\infty(\mathbb{R}^n,\mathbb{R})$, so it holds true working in charts.
To conclude, just apply this proposition to $f=l_g$ for all $g\in G$.