RC circuit analysis with parallel capacitances (using Thevenin equivalent)
I believe that simplifying the circuit to its Thevenin equivalent is
  the correct way to proceed but I'm not sure how to


An alternative way is to work in the frequency domain, use an on-line solver that calculates the transfer function and then use "inverse Laplace" to solve the time domain step function. In the frequency domain we have: -

If we plug the values into this solver we get: -



Then go to this inverse Laplace solver, plug in the formula from above: -



Noting that I've multiplied the denominator by "s" because the input is a step function

And the answer is: -

$$-e^{-11000t}\cosh \left(\sqrt{101000000}t\right)-\frac{11000}{\sqrt{101000000}}e^{-11000t}\sinh \left(\sqrt{101000000}t\right)+1$$

This needs to be multiplied by 9 because the input step voltage is 9 volts and not 1 volt.


  I am trying to determine the voltage across capacitor C3 (voltmeter
  VM1) with respect to time.


If you need a formulaic answer then use the one above. Or, just draw the circuit in your favourite sim tool to get something like this: -



And just in case the formula derived doesn't stack up, here's a double plot of Vout and the derived formula (with the derived formula being made 1 volt higher so the two can be compared): -



Note the "10" in the formula indicated by the arrow is "9+1", just to shift it a tad higher than Vout.