There is a way to derive that formula without using ray-transfer matrices, but instead using the lens equation. First, there is nothing wrong in the way you write $u_2$.

For the first ($f_1$) and second lens ($f_2$), separated by a distance $d$, it holds

$\frac{1}{f_1}=\frac{1}{s_1}+\frac{1}{s_m} \tag{1}$

and

$\frac{1}{f_2}=\frac{1}{s_2}+\frac{1}{d-s_m} \tag{2}$,

where $s_m$ is the position of the image from $s_1$ formed with respect to lens 1. The final image is formed at distance $s_2$ after the second lens.

The piece of information missing in this analysis is that you must leave some distance in front $d_f$ and behind $d_b$ the equivalent lens in order to make things work. The equation for the effective focal length is therefore

$\frac{1}{d_b+s_2}+\frac{1}{d_f+s_1}=\frac{1}{f_e}$

or, rewriting:

$\frac{1}{d_b+s_2}+\frac{1}{d_f+s_1}-\frac{1}{f_e} = 0 \tag{3}$

The calculation proceeds as follows:

1. Write $s_2(s_m)$ as a function of $s_m$ using equation (2) and substitute it in equation (3)
2. Write $s_m(s_1)$ as a function of $s_1$ using equation (1) and substitute it in equation (3). Now eq. (3) features only $s_1,d_b,d_f,f_e$.
3. Write the resulting equation (3) as a fraction. Assuming the denominator is not $0$, we can solve for the numerator $= 0$. This numerator happens to be a quadratic equation in $s_1$ that is $ a_2 s_1^2 + a_1 s_1 + a_0 = 0$ for any value of $s_1$. A quadratic polynomial is always $0$ iff its three coefficients are $ a_2 = a_1 = a_0 = 0$.
4. You need to solve for $f_e, d_f, d_b$ using the three equations given by $ a_2 = a_1 = a_0 = 0$. Now you realize why solving for only one parameter does not work, because the system would be over-constrained.

You will get the two-lens formula for $f=f_e$, the effective focal length, given in the OP.

There is a way to derive that formula without using ray-transfer matrices, but instead using the lens equation. There is nothing wrong in the way you write $u_2$.

For the first ($f_1$) and second lens ($f_2$), separated by a distance $d$, it holds

$\Large \frac{1}{f_1}=\frac{1}{s_1}+\frac{1}{s_m} \tag{1}$

and

$\Large \frac{1}{f_2}=\frac{1}{d-s_m}+\frac{1}{s_2} \tag{2}$,

where $s_m$ is the position of the image from $s_1$ formed with respect to lens 1. The final image is formed at distance $s_2$ after the second lens.

The piece of information missing in this analysis is that you must leave some distance $d_f$ in front and $d_b$ behind the equivalent lens in order to make things work. The equation for the effective focal length is therefore

$\Large \frac{1}{d_f+s_1}+\frac{1}{d_b+s_2}=\frac{1}{f_e}$

or, rewriting:

$\Large \frac{1}{d_f+s_1}-\frac{1}{f_e}+\frac{1}{d_b+s_2} = 0 \tag{3}$

The calculation proceeds as follows:

1. Write $s_2(s_m)$ as a function of $s_m$ using equation (2) and substitute it in equation (3)
2. Write $s_m(s_1)$ as a function of $s_1$ using equation (1) and substitute it in equation (3). Now eq. (3) features only $s_1,d_b,d_f,f_e$.
3. Write the resulting equation (3) as a fraction. Assuming the denominator is not $0$, we can solve for the numerator $= 0$. This numerator happens to be a quadratic equation in $s_1$ that is $ a_2 s_1^2 + a_1 s_1 + a_0 = 0$ for any value of $s_1$. A quadratic polynomial is always $0$ iff its three coefficients are $ a_2 = a_1 = a_0 = 0$.
4. You need to solve for $f_e, d_f, d_b$ using the three equations given by $ a_2 = a_1 = a_0 = 0$. Now you realize why solving for only one parameter does not work, because the system would be over-constrained.

You will get the two-lens formula for $f=f_e$, the effective focal length, given in the OP.

There is a way to derive that formula without using ray-transfer matrices, but instead using the lens equation. First, there is nothing wrong in the way you write $u_2$.

For the first ($f_1$) and second lens ($f_2$), separated by a distance $d$, it holds

$\frac{1}{f_1}=\frac{1}{s_1}+\frac{1}{s_m} \tag{1}$

and

$\frac{1}{f_2}=\frac{1}{s_2}+\frac{1}{d-s_m} \tag{2}$,

where $s_m$ is the position of the image from $s_1$ formed with respect to lens 1. The final image is formed at distance $s_2$ after the second lens.

The piece of information missing in this analysis is that you must leave some distance in front $d_f$ and behind $d_b$ the equivalent lens in order to make things work. The equation for the effective focal length is therefore

$\frac{1}{d_b+s_2}+\frac{1}{d_f+s_1}=\frac{1}{f_e}$

or, rewriting:

$\frac{1}{d_b+s_2}+\frac{1}{d_f+s_1}-\frac{1}{f_e} = 0 \tag{3}$

The calculation proceeds as follows:

1. Write $s_2(s_m)$ as a function of $s_m$ using equation (2) and substitute it in equation (3)
2. Write $s_m(s_1)$ as a function of $s_1$ using equation (1) and substitute it in equation (3). Now eq. (3) features only $s_1,d_b,d_f,f_e$.
3. Write the resulting equation (3) as a fraction. Assuming the denominator is not $0$, we can solve for the numerator $= 0$. This numerator happens to be a quadratic equation in $s_1$ that is $ a_2 s_1^2 + a_1 s_1 + a_0 = 0$ for any value of $s_1$. A quadratic polynomial is always $0$ iff its three coefficients are $ a_2 = a_1 = a_0 = 0$.
4. You need to solve for $f_e, d_f, d_b$ using the three equations given by $ a_2 = a_1 = a_0 = 0$. Now you realize why solving for only one parameter does not work, because the system would be over-constrained.

You will get the two-lens formula for $f_e$ given in the OP.

There is a way to derive that formula without using ray-transfer matrices, but instead using the lens equation. First, there is nothing wrong in the way you write $u_2$.

For the first ($f_1$) and second lens ($f_2$), separated by a distance $d$, it holds

$\frac{1}{f_1}=\frac{1}{s_1}+\frac{1}{s_m} \tag{1}$

and

$\frac{1}{f_2}=\frac{1}{s_2}+\frac{1}{d-s_m} \tag{2}$,

where $s_m$ is the position of the image from $s_1$ formed with respect to lens 1. The final image is formed at distance $s_2$ after the second lens.

The piece of information missing in this analysis is that you must leave some distance in front $d_f$ and behind $d_b$ the equivalent lens in order to make things work. The equation for the effective focal length is therefore

$\frac{1}{d_b+s_2}+\frac{1}{d_f+s_1}=\frac{1}{f_e}$

or, rewriting:

$\frac{1}{d_b+s_2}+\frac{1}{d_f+s_1}-\frac{1}{f_e} = 0 \tag{3}$

The calculation proceeds as follows:

1. Write $s_2(s_m)$ as a function of $s_m$ using equation (2) and substitute it in equation (3)
2. Write $s_m(s_1)$ as a function of $s_1$ using equation (1) and substitute it in equation (3). Now eq. (3) features only $s_1,d_b,d_f,f_e$.
3. Write the resulting equation (3) as a fraction. Assuming the denominator is not $0$, we can solve for the numerator $= 0$. This numerator happens to be a quadratic equation in $s_1$ that is $ a_2 s_1^2 + a_1 s_1 + a_0 = 0$ for any value of $s_1$. A quadratic polynomial is always $0$ iff its three coefficients are $ a_2 = a_1 = a_0 = 0$.
4. You need to solve for $f_e, d_f, d_b$ using the three equations given by $ a_2 = a_1 = a_0 = 0$. Now you realize why solving for only one parameter does not work, because the system would be over-constrained.

You will get the two-lens formula for $f=f_e$, the effective focal length, given in the OP.