Your system has 1 input, the force $F$, and two outputs, the position of $m_1$ and $m_2$. You should end with two equations, one which describes the behavior of $F$ to $m_1$ and one which describes the behavoir from $F$ to $m_2$, i.e. $\frac{U_1(s)}{F(s)}$ and $\frac{U_2(s)}{F(s)}$. These two equations then describe the behavoir of the coupled system.
So your initial laplace transformed equations, I believe, look something like $\alpha U_1(s) + F(s) = \beta U_2(s)$ and $\gamma U_1(s) = \beta U_2(s)$, where $\alpha, \beta, \gamma$ are polynomials of $s$ and the system parameters. Then you just first rewrite the last equations to $U_1(s) = \beta/\gamma U_2(s)$ such that you can subsitute it in the first by which you obtain $U_1(s)/F(s)$. Then rewrite to $U_2(s) = \gamma/\beta U_1(s)$ and subsitute again to obtain $U_2(s)/F(s)$.
You can simulate the behavoir simply with Simulink or Matlab itself e.g. step(tf) or lsim(tf) with tf the transfer function $U_1(s)/F(s)$ or $U_2(s)/F(s)$.
--edit--
Full solution...
We have the following
$$
m_1\ddot{x}_1 = k_2(x_2 - x_1) + d_2(\dot{x}_2 - \dot{x}_1) - k_1x_1 - d_1\dot{x}_1 \\
m_2\ddot{x}_2 = F - k_2(x_2 - x_1) - d_2(\dot{x}_2 - \dot{x}_1)$$
Lets transfer this to the laplace domain and group the terms together such $\alpha X_1(s) = \beta X_2(s)$ and $\gamma X_2(s) = F + \zeta X_1(s)$. With $\alpha, \beta, \gamma, \zeta$ polynomials of $s$ and the coefficients $m_i,k_i,d_i$. We get
$$
\left(m_1s^2 + (d_1 + d_2)s + (k_1 + k_2)\right)X_1(s) = \left(d_2s + k_2\right)X_2(s) \\
\left(m_2s^2 + d_2s + k_2\right)X_2(s) = F(s) + \left(d_2s + k_2\right)X_1(s) \Leftrightarrow F(s) = \left(m_2s^2 + d_2s + k_2\right)X_2(s) - \left(d_2s + k_2\right)X_1(s)
$$
Rewrite the first equation to $\alpha X_1(s) = \beta X_2(s)$ to $X_1(s) = \frac{\beta}{\alpha}X_2(s)$ and $X_2(s) = \frac{\alpha}{\beta}X_1(s)$ we get
$$
X_1(s) = \frac{d_2s + k_2}{m_1s^2 + (d_1 + d_2)s + (k_1 + k_2)}X_2(s) \\
X_2(s) = \frac{m_1s^2 + (d_1 + d_2)s + (k_1 + k_2)}{d_2s + k_2}X_1(s)
$$
Now lets plugin each one such that we retrieve expressions for $\frac{X_1(s)}{F(s)}$ and $\frac{X_2(s)}{F(s)}$. Lets do $\frac{X_1(s)}{F(s)}$ first, we get
$$
F(s) = \left(m_2s^2 + d_2s + k_2\right)\frac{m_1s^2 + (d_1 + d_2)s + (k_1 + k_2)}{d_2s + k_2}X_1(s) - \left(d_2s + k_2\right)X_1(s) \\
F(s) = \frac{\left(m_2s^2 + d_2s + k_2\right)\left(m_1s^2 + (d_1 + d_2)s + (k_1 + k_2)\right) - \left(d_2s + k_2\right)^2}{d_2s + k_2}X_1(s)
$$
Now as you can see the numerator term is a quiet nasty one... You have to first multiply everything and then group all terms such that you get a polynomial of $s$ and the coefficients. For this reason it is also better now to start filling in the values of the coefficients which will make this job then much easier. If you want a analytic solution you can 1) compute it yourself by multiplying everything and collecting the terms 2) let Matlab do this for you. I choose to do it with Matlab.
syms m_1 m_2 d_1 d_2 k_1 k_2 s;
simple((m_2*s^2 + d_2*s + k_2)*(m_1*s^2 + (d_1 + d_2)*s + (k_1 + k_2)) - (d_2*s + k_2)^2)
The simple command will call several commands which could simplify the equation. In this case you will always need the answer for the collect function because this will group the coefficients. In this case it will return
$$
\kappa = (m_1m_2)s^4 + (d_2m_1 + m_2(d_1 + d_2))s^3 + (m_2(k_1 + k_2) + k_2m_1 - d_2^2 + d_2(d_1 + d_2))s^2 + (k_2(d_1 + d_2) - 2d_2k_2 + d_2(k_1 + k_2))s + k_2(k_1 + k_2) - k_2^2
$$
Hence you will get
$$
\frac{X_1(s)}{F(s)} = \frac{d_2s + k_2}{\kappa}
$$
I hope it is clear now and you can solve $\frac{X_2(s)}{F(s)}$ yourself.
-- edit oh well lets do $X_2(s)$ also... --
$$
F(s) = \left(m_2s^2 + d_2s + k_2\right)X_2(s) - \left(d_2s + k_2\right)\frac{d_2s + k_2}{m_1s^2 + (d_1 + d_2)s + (k_1 + k_2)}X_2(s) \\
F(s) = \left(m_2s^2 + d_2s + k_2\right)X_2(s) - \frac{(d_2s + k_2)^2}{m_1s^2 + (d_1 + d_2)s + (k_1 + k_2)}X_2(s)
F(s) = \frac{\left(m_1s^2 + (d_1 + d_2)s + (k_1 + k_2)\right)\left(m_2s^2 + d_2s + k_2\right) - \left(d_2s + k_2\right)^2}{m_1s^2 + (d_1 + d_2)s + (k_1 + k_2)}X_2(s)
\frac{X_2(s)}{F(s)} = \frac{m_1s^2 + (d_1 + d_2)s + (k_1 + k_2)}{\kappa}
$$
Best Answer
My initial guess for modelling the forces involved in firing a rifle would be to model the rifle, bullet and pressurized gas as a piston in which the gas would experience ideal adiabatic expansion until the bullet leaves the barrel. The force exerted by the pressure of the gas on the barrel in the radial direction will cancel out. The pressure does exerts a force on the riffle and the bullet, both proportional to the cross sectional area, $A$, of the barrel (the base of the cylindrical volume of the gas in the barrel). These forces will contribute to the acceleration of both objects, which will result in a increase of the enclosed volume of the gas and therefore a drop in pressure. When assuming the only force involved is due to pressure, the following differential equations can be composed: $$ p=p_0\left(1+\frac{A}{V_0}(x_b-x_r)\right)^{-\gamma}, $$ $$ \frac{\partial^2}{\partial t^2}x_b=\frac{pA}{m_b}=\frac{p_0A}{m_b}\left(1+\frac{A}{V_0}(x_b-x_r)\right)^{-\gamma}, $$ $$ \frac{\partial^2}{\partial t^2}x_r=-\frac{pA}{m_r}=-\frac{p_0A}{m_r}\left(1+\frac{A}{V_0}(x_b-x_r)\right)^{-\gamma}, $$ where $x_b$ and $x_r$ are the positions of the bullet and rifle, $m_b$ and $m_r$ the masses of the bullet and rifle, $p$ is the pressure of the gas, $p_0$ is the starting pressure of the gas, $V_0$ is the starting volume of the pressurized gas and $\gamma$ is the adiabatic index.
This model does has quite some simplifications, such as no friction between barrel and bullet, the air in front of the bullet in the barrel has no influence and no high pressure gas can escape between the bullet and the barrel.
However even this simplified model might be to complex and unnecessary. Because by knowing the muzzle velocity and using the conservation of momentum you that you would have to dissipate the same impulse as how much the bullet has gained, but distribute it over a longer period to reduce the the average force.
Ideally you would want some kind of friction between the moving and static part of the rifle. Because this would mean that there is a constant force acting between the two part and thus on you and therefore minimizing the maximum force experienced and the distance needed to transfer the impulse of the moving part to the user. However friction often means wear, so other systems might be a better option if smallest distance is not that important.
To answer your second problem I would have to simplify it. Say that after firing a bullet of mass $m_b$ at velocity $v_b$, the barrel (and other parts of the rifle moving along with it) with mass $m_r$ should have an equal momentum but in the opposite direction as that of the bullet. From this you can find the initial speed, $v_r$, of the moving part of the rifle: $$ v_r=v_b\frac{m_b}{m_r} $$ After this the spring and the damper will exert a force on this moving part, however the same force will also be exerted onto you: $$ F=m_r\ddot{x}=-kx-c\dot{x} $$ This equation has a known solution and by using the initial conditions $x(t=0)=0$ and $\dot{x}(t=0)=v_r$ you can find the equation of motion and thus the exerted force.