For me, the original data looks to have like a decreasing or constant trend but stl() is giving a different trend altogether. Can someone here please explain why? The decomposition plot is as follows:
To my eye, the second half of 2018 is substantially higher than the prior two years, so the trend that stl() is giving seems not unreasonable.
Trying to fit a time series model to only three years of data is going to require you to take a pretty strong hand guiding the model; I wouldn't expect a turn-key solution to give satisfactory results.
One obvious problem with the decomposition from stl() is that the seasonal figure is wildly overfit. I would try using the Fourier series approach to fitting the seasonal figure that Rob Hyndman describes in this blog post. I have applied that technique in exactly your situation and gotten decent results.
The closest to a turn-key solution is probably going to Facebook's prophet library. But with only three years of data, I still suspect you'll overfit the seasonal component if you call it with the defaults.
If you post your actual data and describe your goals more completely, folks might be able to help more.