How do I check results in Mathematica?

Mathematica is one of the most reliable pieces of software ever written and every version has passed a great number of tests before it is shipped. But no piece of software is And, on the user side, there may always be typos or errors in inputs. Of course, the first line of defense against incorrect results is common sense. For example, one should always ask: Is the order of magnitude plausible? Is the result plausible? 
Additionally, there are some simple tests you can do in Mathematica to check your results without a lot of effort. In more complex cases, more complex tests might be necessary.

Substitiute the result back into the problem

In many cases, you can simply put the result back into the original problem and see if it really solves it.

For example, take the result to this well known problem:

Solve[a (x^2) + b x  +  c == d,  x]

/. is shortcut notation for ReplaceAll, % refers to the last result Mathematica returned. Using FullSimplify, gives {True, True}, showing that both results returned satisfy the equation.

a (x^2) + b x  +  c == d /. % // FullSimplify

The same methode works for the differential equation solver DSolve and other commands that return rules.

DSolve[((y'[x])^2)+2 y[x]==0,y,x]

FullSimplify[((y'[x])^2)+2 y[x]==0/.%]

You can easily substitue multiple results, too.

eqn={y[x]==-z'[x],z[x]==-y'[x]};

sol = DSolve[eqn, {y, z}, x]

FullSimplify[eqn /. sol]

For many other operations, you can simply use the inverse command. For example,

res = Integrate[Sin[x^2] + Cos[x^3], x]

D[res, x] == Sin[x^2] + Cos[x^3] // FullSimplify

Numerical tests

Substitute the result back into the original problem

Let's compare the numerical solution to a differential equation with machine precisions with a numerical solution when using a  higher precision setting.

lowsol = NDSolve[{y'[x] ==  (1/(2  y[x])), y[ .01] ==  .1}, 
          y, {x,  .01, 1}];

highsol = NDSolve[{y'[x] ==  (1/(2  y[x])), y[ .01] ==  .1}, 
          y, {x,  .01, 1}, WorkingPrecision -> 22];

The following plot shows the logarithm of the difference between the numerical solution and the differential equation with the default setting of WorkingPrecision (black) and with WorkingPrecision set to 22 (red).

Plot[
        Evaluate[-Log[
              Abs[{(y'[x] -  (1/(2  y[x])) ) /.  
                    lowsol, (y'[x] -  (1/(2  y[x])) ) /.  
                    highsol}], 10]], {x,  .01, 1}, 
        PlotStyle -> {GrayLevel[0], RGBColor[1, 0, 0]},  
        AxesOrigin -> {0, 0}];

\(Show[%3, \ PlotRange -> {{0, 1}, {0, 0.35}}, 
      PlotLabel -> 
        StyleForm["Relative Error Using Different Precision Settings", "SubSection"]]\)

<< Graphics`Legend`

Show[%7,  Legend->

As you can see, the numerical errors are significantly smaller when using the higher working precision, the rather high first value is an artifact. The trade-off for using higher precision is, of course, increased calculation times.

Random testing

Often, you can get a pretty good indication of whether a result is plausible by plugging in random numbers.

res = Integrate[Sin[Tan[x]], x]

The following command calculated the difference between the differentiated result and the original integrand at 20 random points.

Table[
      Chop[(D[res, x] - Sin[Tan[x]]) /. 
          x -> Random[Complex]], {20}]