What I recommend doing is finding common denominator before plugging in, like this:

## Node A

$$\frac{V_1-V_{ref}}{R_2} + \frac{V_1-V_{x}}{R_1} + \frac{V_1-V_{2}}{R_f} = 0 $$
$$ V_1(R_1R_2 + R_1R_f + R_2R_f) - V_2R_1R_2 - V_x(R_fR_2) + V_{ref}R_fR_1 = 0$$
$$V_x(R_fR_2) = V_1(R_1R_2 + R_1R_f + R_2R_f) - V_2R_1R_2 + V_{ref}R_fR_1$$

## Node B

$$ \frac{V_2-V_{x}}{R_3} + \frac{V_2-V_{1}}{R_f} + \frac{V_2-V_{out}}{R_4} = 0 $$
$$ V_2(R_3R_4 + R_4R_f + R_3R_f) - V_1R_3R_4 - V_x(R_fR_4) + V_oR_fR_3 = 0$$
$$ V_x(R_fR_4) = V_2(R_3R_4 + R_4R_f + R_3R_f) - V_1R_3R_4 + V_oR_fR_3$$

since $$R_4=R_2 \space and \space R_3=R_1$$
$$= V_x(R_fR_2) = V_2(R_1R_2 + R_2R_f + R_1R_f) - V_1R_1R_2 + V_oR_fR_1$$

Then plug it in:
$$V_1(R_1R_2 + R_1R_f + R_2R_f) - V_2R_1R_2 + V_{ref}R_fR_1 = V_2(R_1R_2 + R_2R_f + R_1R_f) - V_1R_1R_2 + V_oR_fR_1$$
$$V_oR_fR_1 = V_1(R_1R_2 + R_1R_f + R_2R_f) + V_1R_1R_2 - V_2(R_1R_2 + R_2R_f + R_1R_f) - V_2R_1R_2 + V_{ref}R_fR_1$$
I'm sure you can figure out the rest.

So here is your schematic without all the weird arrows that you added:

^{simulate this circuit – Schematic created using CircuitLab}

First step is to pick a node and call it zero (ground.) You get to do this with exactly one node of any schematic. If I see a node with LOTS of branches ending in it, I often will pick that one as ground as it reduces the number of terms in some of the expressions. But it really doesn't matter in the end.

Given how you assigned node names, I've selected this ground:

^{simulate this circuit}

Second step is to dump \$R_1\$. It serves no purpose. A current source has infinite impedance and any finite resistance in series with it is entirely irrelevant. The *only* purpose it serves is to change the current source's compliance voltage (the change of which of course is immediately dropped by the change in the resistance.)

So just dump it.

^{simulate this circuit}

Third step is to label the remaining nodes.

I followed your guidance there:

^{simulate this circuit}

Given that we've both arrived at the same number of nodes, I suspect you already have the skills to apply what I did above. So that's good.

Now, as the fourth step, just write out equations for each of the labeled nodes:

$$\begin{align*}
\frac{V_A}{R_3} &= \frac{V_C}{R_3} + I_{E_2} + I_{A_1}\tag{$V_A$}\\\\
\frac{V_B}{R_4} + \frac{V_B}{R_6} + I_{A_1} &= \frac{V_D}{R_4} + \frac{0\:\textrm{V}}{R_6}\tag{$V_B$}\\\\
\frac{V_C}{R_3} + \frac{V_C}{R_5} &= \frac{V_A}{R_3} + \frac{V_D}{R_5} + I_{E_6}\tag{$V_C$}\\\\
\frac{V_D}{R_4} + \frac{V_D}{R_5} + I_{E_2} &= \frac{V_B}{R_4} + \frac{V_C}{R_5}\tag{$V_D$}
\end{align*}$$

This is done by keeping all the outgoing currents on the left side and all of the incoming currents on the right side. Let's discuss this *jonk's* modified nodal method that you will not often see, but is really a whole lot easier to keep track of.

For each node, "imagine" that you are sitting in the middle of the
node (which is like a flat, small, square floor sitting at some
elevation you don't yet know on a long pole from "down there
somewhere.") Water is spilling from above onto your floor. You don't
know from how high, but you can tell it's spilling down onto your
floor. Also, the water that hits your small floor now flows off the
edges and spills down onto other floors below yours. You don't know
where those are at, either. But you know they must be there. (The
"water" is current, the heights of these small floors/nodes is the
voltage of the floor/node.)

So the left side represents the currents spilling away and off the
edges of your floor and onto other nearby floors. And the right side
represents the currents spilling inward and onto your floor from other
nearby floors.

(In keeping with this mental model, you can imagine that current and
voltage sources may be like conveyors that may move water back upward to
higher floors, so to speak. Also, no water is ever lost and all of it
remains in play and the water cannot accumulate on the floors -- the net
flow onto a floor must equal the net flow off of it -- KCL! -- and hence
the justification for setting the left side equal to the right side.)

I've chosen a "direction" for \$I_{E_2}\$ and \$I_{E_6}\$, made evident by whether I placed that current on the left or right side of an equation. The only rule is that I am consistent about it. Here, I've decided that the current flows out of the (+) terminal of a voltage source.

Now, it's at this point that it would appear that there are six unknowns, \$V_A\$, \$V_B\$, \$V_C\$, \$V_D\$, \$I_{E_2}\$, and \$I_{E_6}\$, and only four equations. One way to resolve this is to have used so-called "supernodes" in the above equation development. But that just introduces another bit of terminology I'd like to avoid here. So let's stick with the above equations.

What else do we know? How about the fact that \$V_A=V_D+V_{E_2}\$ and that \$V_C=V_{E_6}\$? Nice. Now we have six unknowns and six equations. So this works out. We can either substitute in these two equations into the other four, reducing the number of variables to four, or we can keep six unknowns and pile on these two added equations. How you approach this is up to you. But let's perform the substitutions since it yields fewer equations and it's easy to do without much chance for misunderstanding:

$$\begin{align*}
\frac{V_D+V_{E_2}}{R_3} &= \frac{V_{E_6}}{R_3} + I_{E_2} + I_{A_1}\tag{$V_A$}\\\\
\frac{V_B}{R_4} + \frac{V_B}{R_6} + I_{A_1} &= \frac{V_D}{R_4}\tag{$V_B$}\\\\
\frac{V_{E_6}}{R_3} + \frac{V_{E_6}}{R_5} &= \frac{V_D+V_{E_2}}{R_3} + \frac{V_D}{R_5} + I_{E_6}\tag{$V_C$}\\\\
\frac{V_D}{R_4} + \frac{V_D}{R_5} + I_{E_2} &= \frac{V_B}{R_4} + \frac{V_{E_6}}{R_5}\tag{$V_D$}
\end{align*}$$

Now the unknowns are just \$V_B\$, \$V_D\$, \$I_{E_2}\$, and \$I_{E_6}\$.

The matrix solution is then:

$$\begin{align*}
\begin{bmatrix}V_B \\ V_D \\ I_{E_2}\\I_{E_6}\end{bmatrix} &=\begin{bmatrix}0 & \frac{1}{R_3} & -1 & 0 \\ \frac{1}{R_4}+\frac{1}{R_6} & \frac{-1}{R_4} & 0 & 0 \\ 0 & \frac{-1}{R_3}+\frac{-1}{R_5} & 0 & -1 \\ \frac{-1}{R_4} & \frac{1}{R_4}+\frac{1}{R_5} & 1 & 0\end{bmatrix}^{-1}
\begin{bmatrix}\frac{V_{E_6}-V_{E_2}}{R_3} + I_{A_1} \\ - I_{A_1} \\ \frac{V_{E_2}-V_{E_6}}{R_3} - \frac{V_{E_6}}{R_5}\\\frac{V_{E_6}}{R_5}\end{bmatrix}
\end{align*}$$

If you stick those equations into some solver, or do it manually, you should get:

$$\begin{align*}
V_B &=-\frac{5}{8}\:\textrm{V}\\
V_D &=8\frac{3}{4}\:\textrm{V}\\
I_{E_2} &=-8\frac{1}{8}\:\textrm{A}\\
I_{E_6} &=-\frac{5}{8}\:\textrm{A}
\end{align*}$$

And from there you should be able to answer any quantitative questions about the rest of the circuit, such as the value for \$V_A-V_B\$.

## Best Answer

I really think nidhin's answer is the best way to go. Use the Laplace transform, calculate \$I_2(s)\$ and then convert it back into a differential equation if you're not interested in the solution.

I just did it for you to show you what you're up against if you try to avoid Laplace transforms:

The Laplace transform gives the solution (I used a CAS for this):

\$I_2\cdot (a\cdot s^3 + b\cdot s^2 + c\cdot s + d) = V_1\cdot (e\cdot s^3 + f\cdot s^2)\$

Yielding the differential equation of the form:

\$a\frac{d^3i_2}{dt^3} + b\frac{d^2i_2}{dt^2} + c\frac{di_2}{dt} + d\cdot i_2 = e\frac{d^3v_1}{dt^3} + f\frac{d^2v_1}{dt^2}\$

Where the coefficients are (I hope I copied them correctly):

\$a = C_1C_2L(R_3R_4 + R_2R_4 + R_1R_4 + R_1R_2)\$

\$b = C_1C_2R_3(R_2R_4 + R_1R_4 + R_1R_2) + (C_1+C_2)L(R_2 + R_3) + L(C_2R_4 + C_1R_1)\$

\$c = C_2R_3(R_4 + R_2) + C_1R_3(R_2 + R_1)\$

\$d = R_3\$

\$e = C_1C_2L(R_2+R_3+R_4)\$

\$f = C_1(C_2R_2R_3+L)\$

I would not advise doing this without transforming to Laplace. While it is possible to solve it without, it can be very hard to see how the equations need to be substituted.