Determine the covariance of geometric Brownian motion without using stochastic calculus.

Question

I've tried this $$Cov(Gt,Gs)=E[G_tG_s]-E[G_t]E[G_s]$$ For the term $E[G_tG_s]$ $$E[G_tG_s]=E[G(t)(G(s)-G(t)+G(t))]$$ $$=E[G(t)G(s-t)]+E[G(t)^2]$$ $$=E[G(t)]E[G(s-t)]+E[G(t)^2]$$ I know that Mean and second moment of $G_t$ $$E[G(t)]=G_0e^{(\mu+\frac{\sigma^2}{2})t}$$ $$E[G_t^2]=G_0^2e^{2t(\mu+\sigma^2)}$$ I just need to calculate $E[G(s-t)]$ $$E[G(s-t)]=G_0E[e^{\mu(s-t)+\sigma W(s-t)}]$$ $$=G_0e^{\mu(s-t)}E[e^{\sigma W(s-t)}]$$ Where $E[e^{\sigma W(s-t)}]$ $$E[e^{\sigma W(s-t)}]=\int_{-\infty}^{\infty}e^{\sigma x}\frac{1}{\sqrt{2\pi (s-t)}}e^{-\frac{x^2}{2(s-t)}}dx$$ $$=\int_{-\infty}^{\infty}e^{\frac{\sigma^2}{2}(s-t)}e^{-\frac{(x-\sigma (s-t))^2}{2(s-t)}}\frac{1}{\sqrt{2\pi (s-t)}}dx$$ $$=e^{\frac{\sigma^2}{2}(s-t)}$$ Then $$E[G(s-t)]=G_0e^{\mu(s-t)}e^{\frac{\sigma^2}{2}(s-t)}$$ $$=G_0e^{(s-t)(\mu+\frac{\sigma^2}{2})}$$ And so the term $$E[G_tG_s]=G_0e^{(\mu+\frac{\sigma^2}{2})t}G_0e^{(s-t)(\mu+\frac{\sigma^2}{2})}+G_0^2e^{2t(\mu+\sigma^2)}$$ $$=G_0^2e^{(\mu+\frac{\sigma^2}{2})s}+G_0^2e^{2t(\mu+\sigma^2)}$$ Finally $$Cov(G_t,G_s)=G_0^2e^{(\mu+\frac{\sigma^2}{2})s}+G_0^2e^{2t(\mu+\sigma^2)}-G_0^2e^{(\mu+\frac{\sigma^2}{2})(t+s)}$$