Computing the Lie bracket on the Lie group $GL(n, \mathbb{R})$

Perhaps this won't satisfy you, but we know that $ad:\mathfrak{g}\to\mathfrak{gl}(\mathfrak{g})$ and $Ad:G\to GL(\mathfrak{g})$ satisfy $d(Ad)_e=ad$, for any Lie group. Here, since we are working with matrices in $\mathbb{R}^{n^2}$, this implies:

If $X,Y\in\mathfrak{gl}_n$, and $g(t)$ is a path of matrices in $GL_n$ such that $g(0)=I_n$ and $g'(0)=X$, then: $$ ad(X)(Y)=\frac{d}{dt}Ad(g(t))(Y)\bigg|_{t=0}=\frac{d}{dt}g(t)Yg(t)^{-1}\bigg|_{t=0}=XY-YX $$ with the last line being by the product rule.

EDIT: Here's a more hands-on approach, following some notes by Wolfgang Ziller (I always get confused during these calculations! I hope this helps someone.)

Let $X,Y\in T_eG$ be tangent vectors at the identity ($G=GL_n(\mathbb{R}))$. Then since $G$ is an open submanifold of $Mat_n(\mathbb{R})$, we can make a canonical identification of all its tangent spaces with $Mat_n(\mathbb{R})$. Once we do this, if $g\in GL_n(\mathbb{R})$ then $(dL_g)_e(X)=gX$, (since matrix multiplication is linear) so the left invariant vector field associated to $X$ is $gX$ at the point $g\in G$. Similarly the one associated to $Y$ is $gY$ at $g$.

Next, notice that if $y_{ij}$ is the coordinate function on matrices picking out the $(i,j)$th entry in the tangent space $T_{h}G$, then $(gX)(y_{ij})=dy_{ij}(gX)=(gX)_{ij}$ because $y_{ij}$ is a linear function. Note that $(gX)_{ij}$ should be thought of as the function defined in a small neighborhood of $h$ given by $g\mapsto (gX)_{ij}$.

Now we can compute the value of the bracket of $X$ and $Y$ at the identity $e$:

$$ [X,Y](y_{ij})=X(Y(y_{ij}))-Y(X_{ij}))=X(g\mapsto(gY)_{ij})-Y(g\mapsto(gX)_{ij}) $$ Now $d(g\mapsto (gY)_{ij})=\sum\limits_{k=1}^ny_{kj}dy_{ik}$, so $$ X(g\mapsto(gY)_{ij})=\sum\limits_{k=1}^ny_{kj}dy_{ik}(X)=\sum\limits_{k}x_{ik}y_{kj}=(XY)_{ij} $$ Similary, $Y(g\mapsto(gX)_{ij})=(YX)_{ij}$, so $[X,Y](y_{ij})=(XY-YX)_{ij}$. This says that the entries of $[X,Y]$ are the same as the entries of $(XY-YX)$, so we're done.

I've now been able to do the calculation myself with a bit of inspiration from freeRmodule's answer:

Using the calculations presented in the question, we calculate the adjoint operator for the Lie group $G := GL(n, \mathbb{R})$. For $g \in G$ and $x \in M_{n,n}(\mathbb{R})=T_{Id}G$ with corresponding left invariant vector field $X$ we have: $$Ad(g)(x) = ((R_{g^{-1}})_{*}X)(Id) = D_gR_{g^{-1}}(\overbrace{X(g)}^{=gx}) = gxg^{-1}$$

By a statement which holds for Lie groups in general we then calculate the Lie bracket of $x$ and $y$ (noting that the Lie group exponential $exp_G$ on $G$ is simply the matrix exponential $Exp$): \begin{align*} [x,y] &= \frac{d}{dt}|_{t=0} Ad(exp_G(tx))(y) = \frac{d}{dt}|_{t=0}Exp(tx)yExp(-tx) \\ &= \left( \frac{d}{dt}|_{t=0}Exp(tx) \right)\cdot y \cdot Id + Id \cdot y \cdot \left( \frac{d}{dt}|_{t=0}Exp(-tx) \right) \\ &= xy -yx \end{align*}