The problem is that just taking the x[0] and x'[0] of your observed data and plugging them into the differential equation will be extremely sensitive to any kind of noise ; not only the 50Hz hum. You are trusting those two samples (rather than the entirety of your observed data) for your estimation, and you will get incorrect results from that, whether the 50Hz noise is present or not.
How to do it right? First, get rid of the 50Hz noise with a notch filter tuned to 50Hz and maybe to a few harmonics.
Then, what you want to do is fit the parameters $A$ and $\phi$ of a $x(n) = A e^{-\delta n} sin(\omega n + \phi)$ model to your data $y(n) = x(n) + noise$. You already know $\delta$ and $\omega$ since they depend on your known system parameter (damping and oscillation frequency). What you really want to know is the best estimate of $A$ and $\phi$ to reconstruct the "theoretical" part (in proper speak: signal subspace projection) of your observation.
Here is how you can estimate $A$ and $\phi$ without lame methods involving looking at peaks and zerocrossing:
You can write your signal model as $x(n) = Re (\alpha z^n) = \frac{1}{2}\alpha z^n - \frac{1}{2} \bar{\alpha} \bar{z}^n$ where $z$ is the known complex pole $z = e^{-\delta + j\omega}$, $\alpha$ is the complex amplitude $\alpha = Ae^{j \phi}$ and $\bar{z}$ denotes the conjugate.
Lay out a $N\times2$ matrix where N is your number of observation, the first column being $\frac{1}{2} [1, z, \ldots, z^N]$ and the second $-\frac{1}{2} [1, \bar{z}, \ldots, \bar{z}^N]$. Do a least square regression of that matrix (which is a basis of your signal subspace) onto the data matrix $[y(1), \ldots, y(N)]$. This will give you LS estimates of $\alpha$ and $\bar{\alpha}$ as the two solutions. Take the modulus and angle of the first to get $A$ and $\phi$. Plug this back into $A e^{-\delta n} sin(\omega n + \phi)$ to get the best fit theoretical model from your data.
# Estimate amplitude and phase of an exponentially damped model
# given observations y and system parameters delta and omega
n = numpy.arange(y.shape[0])
z = numpy.exp(-n * delta + n * omega * 1j)
z_conj = numpy.vstack((0.5 * z, -0.5 * numpy.conj(z))).T
s, _, _, _ = numpy.linalg.lstsq(z_conj, y)
print "Estimated amplitude: ", numpy.abs(s[0])
print "Estimated phase: ", numpy.angle(s[0])
If you want to also estimate $\delta$ and $\omega$ the procedure is a tad more complex and involves an EVD or SVD of the correlation matrix of your data to get a signal subspace basis, one more EVD to get the estimate of $z$, then the procedure sketched above to retrieve $A$ and $\phi$. You will have a robust estimate of all the parameters of the theoretical model from your data. Search for ESPRIT or MUSIC.
Note that this approach can be used to estimate the phase / frequency of a mixture of any number of damped sine waves. So for example, you could use this to retrieve the amplitude and phase of the 50Hz hum itself. Just put two additional columns in the matrix used for the least square procedure, with powers of $e^{j\omega_{noise}}$ and their conjugates. This will give the best fit amplitudes and phases of an exponentially damped sine wave (described by omega and delta) + 50 Hz hum model.
