You should start with the strain energy density $\psi$, then define:
$$
C_{ijkl} = \frac{\partial^2 \psi}{\partial \epsilon_{ij}\partial \epsilon_{kl}}, $$
and then define $$
\sigma_{ij} = C_{ijkl} \epsilon_{kl}
$$
The remainder of my answer will be about explaining why you have to do it that way. Firstly it is physical, there really is energy associated with strain, and if there weren't there would not be any stress. Secondly, it is linear exactly because we are considering the changes in energy due to small strains.
But let's go back to the minor symmetries. We need $C_{ijkl}=C_{jikl}$ because otherwise $\sigma_{ij}\neq\sigma_{ji}$ (and then we get arbitrarily large, hence unphysical, angular velocities for smaller and smaller regions). But the other minor symmetry is not required. If someone handed you a random rank four tensor, let's call it $B_{ijkl}$, and called it an elastic tensor and it didn't have the second minor symmetry, you can define $C_{ijkl}=(B_{ijkl}+B_{ijlk})/2$ $D_{ijkl}=(B_{ijkl}-B_{ijlk})/2$ and then $B=C+D$ but when D is contracted with a symmetric rank two tensor (like the strain tensor) it gives zero. So the part of the rank four tensor without that second minor symmetry simply doesn't contribute, as an actor it does nothing (when you think all the elasticity does is give you stress from strain). So you might as well assume your tensor has both minor symmetries because it acts like it has the second ($B$ and $C$ act the same on symmetric tensors) and it has to have the first.
Did I bring that up to be pedantic? No, I brought it up because the same thing happens if you contract the elasticity tensor with a rank four symmetric combination of the strain tensor. The part of the elasticity tensor without the major symmetry doesn't contribute to the strain energy density. So a random tensor needs the first minor symmetry to be physical. But you might as well assume it has the second minor symmetry since it doesn't affect the stress-strain relationship. And you might as well assume it has the major symmetry because the part that doesn't will not contribute to the strain energy density.
But it is the strain energy density that is physical, and how it changes is what elasticity is. So you aren't really deriving these symmetries as much as saying that only the symmetric ones generate the physical things you want, energy when given strain. And a real derivation should start with strain energy density and strain, and then just define elasticity from that.
Mike Stone is correct. There is no derivation from Newton's laws, and it is just geometry, but I will present it a little differently. Strain angles and rotation angles are how we parameterize all the 3x3 matrices that strain and rotate 3-vectors. Rotations and strains form the group GL(3,R). This is the group of all invertible 3x3 matrices M of real numbers.
We can describe what these transformations do by just talking about the matrices $M$ that are very close to the identity matrix, where all elements in the matrix $\Theta$ are <<1. All these elements are in radians.
$$
M=I+\Theta
$$
$$
\Theta = \begin{bmatrix}
0 & \theta^{12} &-\theta^{13} \\
-\theta^{12} & 0 & \theta^{23} \\
\theta^{13} &-\theta^{23} & 0 \\
\end{bmatrix}_{Asymmetric}
\
+ \begin{bmatrix}
\epsilon^{11} & \epsilon^{12} & \epsilon^{13} \\
\epsilon^{12} & \epsilon^{22} & \epsilon^{23} \\
\epsilon^{13} & \epsilon^{23} & \epsilon^{33} \\
\end{bmatrix}_{Symmetric}
$$
Now apply M to a vector x to get X. We have moved a piece of a body from x to X.
$$
X^i=M^{ij}x^j=(\delta^{ij}+\Theta^{ij})x^j
$$
$$
u^i=(X^i-x^i)=\Theta^{ij} x^j
$$
Where u is the “displacement” of the point. As we move around to different points x in the body, we will get different u ‘s. Differentiating the last equation gives
$$
\frac{\partial u^i}{\partial x^j}=\Theta^{ij}
$$
Therefore we can express the elements of $\Theta$ also in terms the derivatives of the displacements, remembering the strains $\epsilon$ are the symmetric part of $\Theta$ and the rotations $\theta$ are the antisymmetric part of $\Theta$.
$$
\epsilon^{ij}=\frac{1}{2}[\frac{\partial u^i}{\partial x^j}+\frac{\partial u^j}{\partial x^i}]
$$
$$
\theta^{ij}=\frac{1}{2}[\frac{\partial u^i}{\partial x^j}-\frac{\partial u^j}{\partial x^i}]
$$
So far, we have been talking about doing strains and rotations by small angles. Now, by taking products of these small angle transformations, we build all the matrices $M$ in the group for any size elements in $\Theta$. This shows us the matrices $M$ which do macroscopic rotations and strains.
$$
M=\lim_{n\to \infty}(1+\frac{\Theta}{n})^n=e^{\Theta}=I+\Theta+\dfrac{\Theta^2}{2!}+ \dfrac{\Theta^3}{3!}+ \dfrac{\Theta^4}{4!} +…
$$
The $\theta$ are antisymmetric, make M orthogonal ($M^T =M^{-1}$), and leave the length ${x_1}^2 + {x_2}^2+{x_3}^2$ invariant. Because lengths are invariant, the transformations are called rotations. The $\epsilon$ do not leave lengths invariant and are called strains.
The parameter $\theta^{12}$ means rotate the object about axis1 X axis2. That is, put your right thumb perpendicular to the plane formed by axis1 and axis2, such that you fingers would push axis1 into axis 2. Then in this same rotation direction about your thumb rotate the object by $\theta^{12}$ radians.
The parameter $\epsilon^{11}$ means stretch the object by $(1+\epsilon^{11})$ along axis1.
The parameter $\epsilon^{12}$ means parallelepiped the object in the plane containing axis1 and axis2. For example, a square box with its sides initially along axis1 and axis2, becomes a parallelepiped with its sides tilted inward from axis1 and axis2 and its diagonal from the origin stretched. Both sides now make inward angles of $\epsilon^{12}$ radians with their respective axis1 or axis2.
Examples of large rotation and strain matrices are:
Rotate about axis3 by $\theta^{12}$ radians.
$$
M(\theta^{12})=\begin{bmatrix} cos(\theta^{12}) & -sin(\theta^{12}) & 0 \\ sin(\theta^{12}) & cos(\theta^{12}) & 0 \\ 0 & 0 & 1 \end{bmatrix}
$$
Strain (parallelepiped) about the axis3 by $\epsilon^{12}$ radians.
$$
M(\epsilon^{12})=\begin{bmatrix} cosh(\epsilon^{12}) & sinh(\epsilon^{12}) & 0 \\ sinh(\epsilon^{12}) & cosh(\epsilon^{12}) & 0 \\ 0 & 0 & 1 \end{bmatrix}
$$
Strain (stretch) along the axis1 by $\epsilon^{11}$ radians.
$$
M(\epsilon^{11})=\begin{bmatrix} e^{\epsilon^{11} } & 0 \\ 0 & 1 & 0 \\ 0 & 0 & 1 \end{bmatrix}
$$
Notice that for $\epsilon^{11}<<1$, the fractional change of the length of an object in the 1 direction is $\epsilon^{11}$.
The purpose of this rather too long answer was to show the strains and rotations are intimately related and are a group of transformations we can do with our fingers to any piece of material. How the displacements u change as you move around in the body is just what the transformations cause, but is not the fundamental concept of strain or rotation.
Best Answer
Your original math was correct. The strain energy density should have those factors of two in your original answer, when defined in terms of the tensorial definitions of the shear strains. $$ U=\frac{1}{2}\sigma_{ij}\epsilon_{ij}=\frac{1}{2} ( \sigma_{11}\epsilon_{11} + \sigma_{22}\epsilon_{22} + \sigma_{33}\epsilon_{33} + 2\sigma_{12}\epsilon_{12} + 2\sigma_{13} \epsilon_{13}+ 2\sigma_{23}\epsilon_{23} ) $$ The key is to realize that in switching from tensorial notation:
$$ U=\frac{1}{2}\sigma_{ij}\epsilon_{ij}=\frac{1}{2} \left[ \begin{array}{ccc} \sigma_{11} & \sigma_{12} & \sigma_{13} \\ \sigma_{21} & \sigma_{22} & \sigma_{23}\\ \sigma_{31} & \sigma_{32} & \sigma_{33} \end{array} \right]\ \left[ \begin{array}{ccc} \epsilon_{11} & \epsilon_{12} & \epsilon_{13} \\ \epsilon_{21} & \epsilon_{22} & \epsilon_{23}\\ \epsilon_{31} & \epsilon_{32} & \epsilon_{33} \end{array} \right]\ $$
to engineering (i.e. Voigt) notation, one must account for a change in definition of the shear strains. Explicitly, in engineering notation, we write: $$ \left[ \begin{array} {ccc} \sigma_{1} \\ \sigma_{2} \\ \sigma_{3} \\ \sigma_{4} \\ \sigma_{5} \\ \sigma_{6} \end{array} \right]\ = C_{ij} \left[ \begin{array} {ccc} \epsilon_{1} \\ \epsilon_{2} \\ \epsilon_{3} \\ \epsilon_{4} \\ \epsilon_{5} \\ \epsilon_{6} \end{array} \right]\, where \left[ \begin{array} {ccc} \epsilon_{1} \\ \epsilon_{2} \\ \epsilon_{3} \\ \epsilon_{4} \\ \epsilon_{5} \\ \epsilon_{6} \end{array} \right]\ = \left[ \begin{array} {ccc} \epsilon_{11} \\ \epsilon_{22} \\ \epsilon_{33} \\ 2\epsilon_{23} \\ 2\epsilon_{13} \\ 2\epsilon_{12} \end{array} \right]\ = \left[ \begin{array} {ccc} \epsilon_{11} \\ \epsilon_{22} \\ \epsilon_{33} \\ \gamma_{23} \\ \gamma_{13} \\ \gamma_{12} \end{array} \right]\ $$
In engineering notation, we can also write: $U=\frac{1}{2}\sigma_i \epsilon_i$, and this is equivalent to the tensorial definition; you get the same factors of two when you switch from $\epsilon_i$ to $\epsilon_{ij}$ notation. The various defintions of shear strains can get confusing!