If we assume a perfect two-body problem, absent perturbations from external bodies or non-spherical gravity sources (i.e., perfect conic orbits with no precession or variation), your constraints regarding inclination are actually unnecessary, as we may, without loss of generality, examine this problem in a perifocal reference frame.
Perifocal coordinates
Let the origin of the system be the focus of the orbit, and let the orthogonal unit vectors $\mathbf{\hat{p}}$ and $\mathbf{\hat{q}}$ lie in the orbital plane, such that the periapsis of the orbit is located on a line extending from the origin in the $\mathbf{\hat{p}}$ direction. The unit vector $\mathbf{\hat{w}} = \mathbf{\hat{p}} \times \mathbf{\hat{q}}$ completes the coordinate system. (The orbital angular momentum vector lies in the $\mathbf{\hat{w}}$ direction).
If we further constrain ourselves to only look at periapsis and -- provided $e < 1$ -- apoapsis, our position vector is then strictly in the $\pm \mathbf{\hat{p}}$ direction, and our velocity vector is strictly in the $\pm \mathbf{\hat{q}}$ direction.
Keplerian elements in relation to perifocal coordinates
The Keplerian elements $i$ (inclination), $\Omega$ (argument of ascending node), and $\omega$ (argument of periapsis) are simply Euler angles used to translate between a perifocal coordinate system and an equatorial coordinate system. By focusing on a perifocal coordinate system, we can temporarily dispense with these, leaving only $e$ (eccentricity), $a$ (semi-major axis), and $\nu$ (true anomaly).
With that established, we need to develop two more scalar terms that remain constant in an unperturbed orbit (note: $\mu$ is the gravitational parameter of the central body):
- Specific energy: $$\varepsilon = \frac{v^2}{2} - \frac{\mu}{\|\mathbf{r}\|}$$
- Specific angular momentum: $$h = \|\mathbf{r} \times \mathbf{v}\|$$
To facilitate the upcoming derivations, we will expand these in terms of components (including ones known to be zero):
$$\varepsilon = \frac{v^2_p + v^2_q + v^2_w}{2} - \frac{\mu}{\sqrt{r^2_p + r^2_q + r^2_w}}$$
$$h = \sqrt{(r_q v_w - r_w v_q)^2 + (r_w v_p - r_p v_w)^2 + (r_p v_q - r_q v_p)^2}$$
We now note that
$$a = -\frac{\mu}{2\varepsilon}$$
and $$e = \sqrt{1 + \frac{2 \varepsilon h^2}{\mu^2}}$$
Finally, to your question
To frame this explicitly, you are looking for how the elements are affected by velocity perturbations. This means we want some expressions for these as a function of velocity, all other terms (in particular, the position) being constant.
In general, you want to accomplish this through a multivariate Taylor expansion.
Let $Q : \mathbb{R}^3 \rightarrow \mathbb{R}$ be an arbitrary real-valued continuous and at-least-twice differentiable (since you want second order accuracy) function of its arguments. Then to second order accuracy, you have
$$ Q(\mathbf{x}) = Q(\mathbf{a}) + \Big[(\mathbf{x}-\mathbf{a}) \cdot \mathbf{\nabla}Q(\mathbf{a})\Big] + \frac{1}{2}\Big[(\mathbf{x}-\mathbf{a}) \cdot \mathbf{H}(\mathbf{a}) \cdot (\mathbf{x}-\mathbf{a})\Big] + \epsilon(\|\mathbf{x}-\mathbf{a}\|^3)$$
where $\mathbf{\nabla}Q$ is the gradient of $Q$, i.e., $\mathbf{\nabla}Q = [Q_1 \; Q_2 \; Q_3]^T$, (using subscript numerals to denote partial differentiation with respect to that argument), and $\mathbf{H}$ is the Hessian matrix:
$$H = \begin{bmatrix}Q_{11} & Q_{12} & Q_{13}\\ Q_{21} & Q_{22} & Q_{23}\\ Q_{31} & Q_{32} & Q_{33} \end{bmatrix}$$
Writing this explicitly gets cumbersome quickly and would not be easily readable, especially since it will involve multiple invocations of the chain rule at higher orders in multiple dimensions.
Instead, we will demonstrate one specific case: the semimajor axis (one of the simplest elements to calculate), $a$, to first-order accuracy only.
Let a prefix $\delta$ denote an infinitesimal perturbation to a quantity.
We now have,
$$\delta a = \frac{\partial a}{\partial \varepsilon}\left(\frac{\partial \varepsilon}{\partial v_p}\delta v_p + \frac{\partial \varepsilon}{\partial v_q}\delta v_q + \frac{\partial \varepsilon}{\partial v_w} \delta v_w\right) + \epsilon(\|\delta \mathbf{v}\|^2)$$
Expanding this, we have
$$\delta a = \frac{\mu}{2\varepsilon^2} \Big(v_p \delta v_p + v_q \delta v_q + v_w \delta v_w \Big) + \epsilon(\|\delta \mathbf{v}\|^2)$$
What can we draw from this? Note that the parenthesized term is identically equal to $\mathbf{v} \cdot \mathbf{\delta v}$. That is to say, to first-order accuracy, the semimajor axis is affected by the component of the perturbation in the prograde/retrograde direction only. This is the case regardless of where in the orbit this happens, not just at one of the apsides.
If we want to explore a second order expansion, it gets ugly quickly.
$$\delta a = \frac{\partial a}{\partial \varepsilon}\left(\frac{\partial \varepsilon}{\partial v_p}\delta v_p + \frac{\partial \varepsilon}{\partial v_q}\delta v_q + \frac{\partial \varepsilon}{\partial v_w} \delta v_w\right) + \frac{1}{2}\left\{ \frac{\partial a}{\partial \varepsilon} \left(
\frac{\partial^2 \varepsilon}{\partial v^2_p}(\delta v_p)^2 +
\frac{\partial^2 \varepsilon}{\partial v_p \partial v_q}(\delta v_p \delta v_q) +
\frac{\partial^2 \varepsilon}{\partial v_p \partial v_w}(\delta v_p \delta v_w) +
\frac{\partial^2 \varepsilon}{\partial v_q \partial v_p}(\delta v_q \delta v_p) +
\frac{\partial^2 \varepsilon}{\partial v^2_q}(\delta v_q)^2 +
\frac{\partial^2 \varepsilon}{\partial v_q \partial v_w}(\delta v_q \delta v_w) +
\frac{\partial^2 \varepsilon}{\partial v_w \partial v_p}(\delta v_w \delta v_p) +
\frac{\partial^2 \varepsilon}{\partial v_w \partial v_q}(\delta v_w \delta v_q) +
\frac{\partial^2 \varepsilon}{\partial v^2_w}(\delta v_w)^2
\right)
+ \frac{\partial^2 a}{\partial \varepsilon^2} \left[
\left(\frac{\partial \varepsilon}{\partial v_p}\right)^2(\delta v_p)^2 +
\left(\frac{\partial \varepsilon}{\partial v_p}\right)\left(\frac{\partial \varepsilon}{\partial v_q}\right)(\delta v_p \delta v_q) +
\left(\frac{\partial \varepsilon}{\partial v_p}\right)\left(\frac{\partial \varepsilon}{\partial v_w}\right)(\delta v_p \delta v_w) +
\left(\frac{\partial \varepsilon}{\partial v_q}\right)\left(\frac{\partial \varepsilon}{\partial v_p}\right)(\delta v_q \delta v_p) +
\left(\frac{\partial \varepsilon}{\partial v_q}\right)^2(\delta v_q)^2 +
\left(\frac{\partial \varepsilon}{\partial v_q}\right)\left(\frac{\partial \varepsilon}{\partial v_w}\right)(\delta v_q \delta v_w) +
\left(\frac{\partial \varepsilon}{\partial v_w}\right)\left(\frac{\partial \varepsilon}{\partial v_p}\right)(\delta v_w \delta v_p) +
\left(\frac{\partial \varepsilon}{\partial v_w}\right)\left(\frac{\partial \varepsilon}{\partial v_q}\right)(\delta v_w \delta v_q) +
\left(\frac{\partial \varepsilon}{\partial v_w}\right)^2(\delta v_w)^2
\right]\right\} + \epsilon(\|\delta \mathbf{v}\|^3)$$
Hoo boy!
We now make the following substitutions:
$$\frac{\partial a}{\partial \varepsilon} = \frac{\mu}{2\varepsilon^2}$$
$$\frac{\partial^2 a}{\partial \varepsilon^2} = -\frac{\mu}{4\varepsilon^3}$$
$$\frac{\partial \varepsilon}{\partial v_i} = v_i,\;i=p,q,w$$
$$\frac{\partial^2 \varepsilon}{\partial v_i \partial v_j} = \delta_{ij} = \left\{\begin{array}{l}1, & i=j\\0, & i \neq j\end{array}\right., \; i,j = p,q,w$$
Now, finally
$$\delta a = \frac{\mu}{2\varepsilon^2} \Big(v_p \delta v_p + v_q \delta v_q + v_w \delta v_w \Big) +
+ \frac{1}{2}\left\{ \frac{\mu}{2 \varepsilon^2} \Big[
(\delta v_p)^2 +
(\delta v_q)^2 +
(\delta v_w)^2
\Big]
- \frac{\mu}{4 \varepsilon^3} \Big[
v_p^2(\delta v_p)^2 +
v_q^2(\delta v_q)^2 +
v_w^2(\delta v_w)^2 +
2 v_p v_q(\delta v_p \delta v_q) +
2 v_p v_w(\delta v_p \delta v_w) +
2 v_q v_w(\delta v_q \delta v_w)
\Big]\right\}
+\epsilon(\|\delta \mathbf{v}\|^3)$$
Writing this more compactly:
$$\delta a = \frac{\mu}{2\varepsilon^2} \left( \mathbf{v} \cdot \delta \mathbf{v} + \frac{\|\delta \mathbf{v}\|^2}{2} \right) - \frac{\mu}{4 \varepsilon^3} \Big( \mathbf{v} \cdot \delta \mathbf{v} \Big)^2 +\epsilon(\|\delta \mathbf{v}\|^3)$$
Important consequences: at second-order accuracy, the previous assertion regarding prograde/retrograde is not applicable.
Final thoughts
Though there are a few hiccups involved in applying this to the three orientation elements due to some degeneracy in the rotation formalism defined by Euler angles, this answer should provide all the tools you need to calculate what you are asking. It is a simple -- though tedious -- matter of simply following the Taylor expansion for each of the remaining elements and applying coordinate transformations where appropriate.
This should also provide a healthy illustration for why nobody really does second-order expansions. The equations get nasty very quickly, with very little benefit to offer from a numerical analysis standpoint.
