I) Not all equations of motion (eom) are variational. A famous example is the self-dual five-form in type IIB superstring theory. In classical point mechanics, frictional forces typically lead to non-variational problems.
II) Consider for instance $n$ variable $q^i$ and $n$ eoms,
$$\tag{1} E_i~\approx~ 0, \qquad i~\in~\{1, \ldots, n\}. $$
A simplified version of OP's problem (v3) is the following:
Does there exist an action
$$\tag{2} S[q] ~=~\int{\rm d}t~L$$
such that the Euler-Lagrange derivatives
$$\tag{3} \frac{\delta S}{\delta q^i}~=~E_i $$
precisely become the given $E_i$-functions?
The above restricted problem is relatively easy to answer once and for all, because one may differentiate the known $E_i$-functions to arrive at a set of consistency conditions. Let us for simplicity assume that the functions $E_i=E_i(q)$ do not involve generalized velocities $\dot{q}^i$, accelerations $\ddot{q}^i$, and so forth. Then we may assume that the Lagrangian $L$ does not depend on time derivatives of $q^i$ as well. So the question becomes if
$$\tag{4} \frac{\partial L}{\partial q^i}~=~E_i ? $$
We can collect the information of the eoms in a one-form
$$\tag{5} E~:=~E_i ~{\rm d}q^i.$$
The question rewrites as
$$\tag{6} {\rm d}L~=~E? $$
Hence the Lagrangian $L$ exists if $E$ is an exact one-form.
III) However, the above discussion is in many ways oversimplified. The eoms (1) do not have a unique form! E.g. one may multiply the given $E_i$-functions with an invertible $q$-dependent matrix $A^i{}_j$ such that the eoms (1) equivalently read
$$\tag{7} \sum_{i=1}^n E_i A^i{}_j~\approx~ 0. $$
Or perhaps the system variables $q^i$ should be viewed as a subsystem of a larger system with more dynamical or auxiliary variables?
Ultimately, the main question is whether the eoms have an action principle or not; the particular form of the eoms (that the Euler-Lagrange equations spit out) is not important in this context.
This opens up a lot of possibilities, and it can be very difficult to systematically find an action principle; or conversely, to prove a no-go theorem that a given set of eoms is not variational.
Well, if you have a term like $\partial_\mu \mathcal{J}^\mu$, the divergence theorem lets you convert it into a surface term upon integrating to find the action, and since variations are assumed to vanish at the boundary, this term goes away. The Euler-Lagrange equations don't change because they come from setting the variation of the action to zero.
Example: Suppose you have a lagrangian $\mathcal{L}_0$ and add a divergence to get $\mathcal{L} = \mathcal{L}_0 + \partial_\mu \mathcal{J}^\mu$. Recall that the action is (in your favorite number of dimensions):
$$S = \int dx\ \mathcal{L} = S_0 + \int dx\ \partial_\mu \mathcal{J}^\mu = S_0 + \int dS\ n_\mu \mathcal{J}^\mu$$
Here S_0 is the integral of $\mathcal{L}_0$, and $n^\mu$ the normal vector to your boundary.
The equations of motion are the condition that $\delta S = 0$ to first order whenever we make a variation in $\mathcal{L}$. So:
$$\delta S = \delta S_0 + \int dS\ n_\mu \delta(\mathcal{J}^\mu)$$
But $\mathcal{J}^\mu$ is constructed out of the fields for which you want the equations of motion. Since by hypothesis the variation of the fields at the boundary is zero, so is the variation of $\mathcal{J}^\mu$. The last term vanishes, and we get $\delta S = \delta S_0$. This implies that the e.o.m. don't change, since $\delta S = 0$ is equivalent to $\delta S_0 = 0$.
Best Answer
In some sense, yes it is, and in others it is certainly not.
The fact that in the Einstein field equations: $G_{\mu \nu} = 8 \pi T_{\mu \nu}$, you have complete 'freedom' to define $T_{\mu \nu}$ however you would like (within some kind of exceptions), means you can allow for highly non-trivial curvature dependent terms to be coupled to the stress-energy. See this really interesting paper by Capozziello et al (http://arxiv.org/abs/grqc/0703067) that shows how you can 'bunch' the additional curvature terms for an $f(R)$ gravity into the 'effective' stress-energy tensor. In this way, the uniqueness of what you may call the 'Einstein' field equations does not hold.
However, once one moves to the vacuum the answer is suddently a definitive yes. The Einstein equations are generated by having an action that contains a geometric invariant quantity. If it did not then the principles of general covariance would no longer hold since arbitrary coordinate transformations would necessarily destroy the gauge symmetries. Now, that being said, one can consider the most general possible geometric invariant Lagrangian in this sense, it would look something like.. $\mathcal{L} = f(g_{\mu \nu} R^{\mu \nu}, R^{\mu \nu} R_{\mu \nu}, R^{\alpha \beta \gamma \delta} R_{\alpha \beta \gamma \delta}, C^{\alpha \beta \gamma \delta} C_{\alpha \beta \gamma \delta}, g^{\mu \nu} R^{\alpha \beta} R_{\alpha \beta \mu \nu}...)$ and the list of possibilities goes on. Since you want to replicate the Einstein equations in themselves you can immediatly throw away all but the linear terms in the above. Any of the other variants when you perform the calculus of variations gives you higher than second order derivatives.
Now, keeping only the linear terms you see that the action will necessarily be the Einstein-Hilbert one plus (maybe) some other curvature invariants. Since these will not vanish once you perform the variation, unless they either vanish identically or vanish on the boundary as a divergence term, your field equations will necessarily still be the Einstein-Hilbert ones provided you only keep $\mathcal{L} = R^{\sigma}_{\sigma}$.
A rigorous proof of this would be interesting. I am unaware of any in the literature.