The mistakes you are doing are:
You didn't considered the flux coming from them in between them. You have to take all the flux in all directions coming from them. You should take the gaussian across the surface of the plane otherwise you will get wrong result.
$|\vec E_+|=|\vec E_-|=\frac{\sigma}{2\epsilon_0}$ and not $\frac{-\sigma}{2\epsilon_0}$ for $|\vec E_-|.\space$ $\sigma$ is the magnitude of the charge density.
You are incorrectly adding the fields which gave you $0$ inside. The magnitudes have to be added when directions are same and subtracted when directions are opposite.
This is what we get from Gauss's law: $$\vec{E}=\frac{\sigma}{2\epsilon_0}\hat r$$
where, $$|\vec{E}|=\frac{\sigma}{2\epsilon_0}$$where $\sigma$ is the magnitude of surface charge density
So, outside, if direction of $\vec{E_+}$ is $\hat r$ then, direction of $\vec{E_-}$ is $-\hat r$ $$\vec{E_+}=\frac{\sigma}{2\epsilon_0}\hat r$$$$\vec{E_-}=\frac{\sigma}{2\epsilon_0}(-\hat r)$$$$\vec{E_+}+\vec{E_-}=\frac{\sigma}{2\epsilon_0}\hat r+\frac{-\sigma}{2\epsilon_0}\hat r$$ $$=0$$
Inside, both $\vec E_+$ and $\vec E_-$ has same direction $\hat r$ $$\vec{E_+}=\frac{\sigma}{2\epsilon_0}\hat r$$$$\vec{E_-}=\frac{\sigma}{2\epsilon_0}\hat r$$$$\vec{E_+}+\vec{E_-}=\frac{\sigma}{2\epsilon_0}\hat r+\frac{\sigma}{2\epsilon_0}\hat r=\frac{\sigma}{\epsilon_0}\hat r$$
Talking in magnitudes, inside, the magnitudes have to be added, $$|\vec E_+|+|\vec E_-|=\frac{\sigma}{2\epsilon_0}+\frac{\sigma}{2\epsilon_0}=\frac{\sigma}{\epsilon_0}$$
outside, they have to be subtracted, $$|\vec E_+|-|\vec E_-|=\frac{\sigma}{2\epsilon_0}-\frac{\sigma}{2\epsilon_0}=0$$
This is an extremely common mistake in introductory EM - from students who actually spend time thinking about the problem, anyway ;-) Use Gauss's law in both cases:
In the case of infinite plates, you do not have the result you give first. A Gaussian cylinder has two disks on either side of the plate, so
$$E_1(2A)=\frac{\sigma A}{\epsilon_0}\rightarrow E_1=\frac{\sigma}{2\epsilon_0}$$
And from superposition you get the total electric field
$$E=\frac{\sigma}{\epsilon_0}$$
You second case is correct, but the charge enclosed by your surface is $Q/2$ relative to the first case (conservation of charge, if you want the same answer you better have the same total charge on the plates), so
$$E_1A=\frac{(\sigma/2) A}{\epsilon_0}\rightarrow E_1=\frac{\sigma}{2\epsilon_0}$$
Which again gets you the same answer when you apply superposition.
Best Answer
Electrostatic refers to the case where the fields are not time dependent. In that case the Maxwell's equations reduce to:
$$\nabla \cdot E =\frac{\rho}{\epsilon_o} \\ \nabla \times E = 0 \implies E=-\nabla \phi \\ \text{then,} \nabla \cdot \nabla \phi = \nabla^2 \phi = -\frac{\rho}{\epsilon_o} $$
The solution to the last equation is:
$$ \phi = \frac{1}{4\pi\epsilon_o} \int \frac{\rho(\mathbf{r'})}{|\mathbf{r}-\mathbf{r}'|}d^3r' $$
which gives the electrostatic electric field equation you have written for point charges.
But if the case is time dependent then you have the Maxwell equation:
$$\nabla \times E = -\frac{\partial B}{\partial t}$$
and you can no longer define $\phi$ as we did above. In this case you have to work with retarded potentials:
$$\psi(\mathbf{r}, t) = \frac{1}{4\pi\epsilon_o} \int \frac{\rho(\mathbf{r'}, t')}{|\mathbf{r}-\mathbf{r}'|}d^3r'$$
where $\psi = \phi, A_x,A_y,A_z$ and $t'=\frac{|\mathbf{r}-\mathbf{r'}|}{c}$. (Note that when $\psi = \phi$ then $\rho$ is the charge density and when $\psi = A_i$ then $\rho = J_i$, i.e., the $i^{th}$ component of current density). Then the fields are given as:
$$ \mathbf{ E} = -\nabla \phi -\frac{\partial \mathbf{A}}{\partial t}\\ \mathbf{B} = \nabla \times \mathbf{A} $$
These time dependent fields turn out to be very complicated expressions and are named the Jefimenko equations.