I've taken a look on these equations.
I know you have a 1D element, but rewriting the equation would eliminate confusion if you want to write elements in 2D and 3D.
Therefore your stiffness matrix should be
\begin{equation} \mathbf{K}{el} = \sum{i}^{n_{gauss}} \mathbf{B}^T \mathbf{E} \mathbf{B} w_{\xi,i} det(J). \end{equation}
for the numerical integration with Gauss. Your stiffness matrix formula used $J$ instead of $det(J)$, which is not wrong. The determinant of the Jacobian is only significant for multidimensional cases.
The formula for your Jacobian should be \begin{equation} J = \sum_{i}^{n} \frac{\partial N_{i}}{\partial \xi}x_{i}. \end{equation}
The Jacobian in the one-dimensional case is defined as $J = \left[ \frac{\partial x}{\partial \xi} \right]$ for \begin{equation} \begin{bmatrix} \frac{\partial N}{\partial \xi} \end{bmatrix} = J \begin{bmatrix} \frac{\partial N}{\partial x} \end{bmatrix} \end{equation}
\begin{equation} \begin{bmatrix} \frac{\partial N}{\partial \xi} \end{bmatrix} = \begin{bmatrix} \frac{\partial x}{\partial \xi} \end{bmatrix} \begin{bmatrix} \frac{\partial N}{\partial x} \end{bmatrix} \end{equation}
The Jacobian in the isoparametric concept approximates the transformation by the derivatives of the shape function. 526
Here $\frac{\partial x}{\partial \xi}$ has to be approximated. Like the x-coordinate in the element would be calculated by \begin{equation} x = \begin{bmatrix} N_1 & N_2 & N_3 & N_4 \end{bmatrix} \begin{bmatrix} \hat{x}_1 \ \hat{x}_2 \ \hat{x}_3 \ \hat{x}_4 \ \end{bmatrix} \end{equation} $\frac{\partial x}{\partial \xi}$ is approximated by \begin{equation} \frac{\partial x}{\partial \xi} = \begin{bmatrix} \frac{\partial N_1}{\partial \xi} & \frac{\partial N_2}{\partial \xi} & \frac{\partial N_3}{\partial \xi} & \frac{\partial N_4}{\partial \xi} \end{bmatrix} \begin{bmatrix} \hat{x}_1 \ \hat{x}_2 \ \hat{x}_3 \ \hat{x}_4 \ \end{bmatrix}. \end{equation}
Important: The definition of the Jacobian in many books covering finite elements is not consistent to the mathematical definition. So be careful when you are doing your transformations.
When you are doing the transformation of the derivatives of the shape function, you should have \begin{equation} \begin{bmatrix} \frac{\partial N_1}{\partial x_1} & \frac{\partial N_2}{\partial x_2} & \frac{\partial N_3}{\partial x_3} & \frac{\partial N_4}{\partial x_4} \end{bmatrix} = \mathbf{J}^{-1} \begin{bmatrix} \frac{\partial N_1}{\partial \xi_1} & \frac{\partial N_2}{\partial \xi_2} & \frac{\partial N_3}{\partial \xi_3} & \frac{\partial N_4}{\partial \xi_4} \end{bmatrix}, \end{equation}
so that seems to be ok, except that you have the matrix transposed.
Your stiffness matrix then is wrong. Note that you have a dyadic product, where in this case you get a $4 \times 4$ matrix out of $1 \times 4$ vectors. 804