Matlab can only deal with transfer functions that are rational ratios of polynomials in $s$ (or $z$, if you're in that domain). So you must approximate. You can either use the responses that @Matt L. has shown you, and convolve by the impulse response at each step, or you can write the partial differential equation (not an ODE) and solve it at each time step, or you can approximate the transfer function as a rational fraction of polynomials in $s$.
For the rest of this answer I'm going to assume that you really mean $G\left ( s \right)$ is some lazy shorthand for a real-world dimensionally sensible form. If so, then you must really mean that $G \left (s \right ) = k e^{-\sqrt{\tau s}}$ where $k$ is in whatevers over whatever elses (what are your inputs and outputs?) and $\tau$ is in seconds so that $\tau s$ can be in Napirs.
I would do this by finding a ratio of polynomials in $s$ that approximates $e^{-\sqrt{\tau s}}$. Having dinged you for not using dimensions, I'd still probably start with $e^{-\sqrt{\theta}}$ and find an approximate polynomial in $\theta$ (use the letter that's most sensible to you), then scale it into a polynomial in $s$.
You want something that has a good-enough fit for you -- and this changes depending on the problem you want to solve. If you're designing a closed-loop controller, a good enough fit is probably one that doesn't deviate by more than five degrees or so in phase and 1dB or less in amplitude over the frequency range at which your closed-loop gain is higher than -10dB or so -- which means that depending on how aggressively you're designing your controller, you need more or less terms in your approximation.
You will, in essence, be designing a linear time invariant ODE that simulates the partial differential equation whose transfer function is $e^{-\sqrt{\tau s}}$.
Whatever else I did, I'd take the impulse and step response of my approximation and compare them to the ones that @Matt L. has shown you, and I'd make sure that they were close enough to be intuitively satisfying.