Result with Thevenin method doesn't equal node method
I have no idea how you got to the second to last equation of yours. However, let's assume you have the circuit with a Thevenin equivalent circuit substituted in only for R2 and I (your \$R_{th}\$ and \$V_{th}\$ are correct).

Several ways to solve this new circuit:

Nodal Analysis






simulate this circuit &ndash; Schematic created using CircuitLab

Assuming current flows "clockwise" (in the direction of the drawn arrow), performing nodal analysis at node e gives:

\begin{gather}
I_{R1} = I_{R2}\\
\frac{V - e}{R_1} = \frac{e - V_{th}}{R_2}\\
R_2 V + R_1 V_{th} = R_1 e + R_2 e\\
e = \frac{R_2 V + R_1 V_{th}}{R_1 + R_2} = \frac{R_2 V + R_1 R_2 I}{R_1 + R_2}
\end{gather}

Which is the same answer as you got without using the Thevenin equivalent.

Superposition

See Weezveez's answer. It's the same result.

Mesh Analysis

\begin{gather}
V - V_{R1} - V_{R2} - V_{th} = 0\\
V - I_1 R_1 - I_1 R_2 - V_{th} = 0\\
I_1 = \frac{V - V_{th}}{R_1 + R_2}\\
e = V_{th} + V_{R2} = V_{th} + I_1 R_2\\
e = I R_2 + \frac{V - V_{th}}{R_1 + R_2} R_2\\
e = I R_2 + \frac{V - I R_2}{R_1 + R_2} R_2\\
e = \frac{R_2 V + R_1 R_2 I}{R_1 + R_2}
\end{gather}

As you can see, all of these methods give the same answer.