Different software packages use different ways to calculate many statistics. Skewness is a good example (see the attached) where many of the packages use different skewness statistics.  There are also different algorithms for the median (e.g., what to do when the median is between two values). The means should be the same across the main packages, but if you are running means for a large number of variables they may deal with missing values differently.

It is often informative to try a few packages. Seeing the differences and then thinking why there are differences is sometimes useful. For complex analyses there are usually differences.

The by-hand methods are often simpler. Often the packages do something to make the estimate a little better, but they sometimes have options for using the "by hand" equation.

If you show your data readers can show you what they find.

Article Problematic standard errors and confidence intervals for ske...

Simple statistics such as mean or median should not differ, no matter what system is used to derive them. There was probably an error somewhere. Most likely in the hand calculations, but possibly in the data entry stage.

Some high-level statistics will differ depending on how they were derived - for example, whether a multiple regression formula was arrived at iteratively or by rotating a matrix. And sometimes a statistical program will be set up to take a "short cut" in computing some parameters - or it may just be programmed incorrectly. But the mean? The median? No way.

I agree with Stephen that the problem is likely with the hand calculation, but I would be less certain about all packages calculating things that seem simple, like the median, identically.

If we define the median is the estimated 50th percentile and use the quantile function in R, it allows 9 ways to calculated these. For example, the median for 1,2,3,4 is 2.5 in R's default but 2 with some of the alternative methods. It is often useful to be able to check how the packages you are using calculate different statistics.  Of course if it is defined explicitly, like the value the minimizes sum(abs(est median - x_i)), then it should produce the same thing.

quantile(1:4,.5,type=7)

50%

2.5

> quantile(1:4,.5,type=3)

50%

2

So in R if you type median.default and look towards the bottom it shows that if there are an even number of cases it takes the mean for the cases just each side of the halfway point.