Numerical instability in cosh and sinh - integral functions

You are calculating the small difference between numbers which are getting quite large and the default WorkingPrecision of Plot (which is MachinePrecision, usually about 16) is just not high enough.

So simply increase the WorkingPrecision of Plot, e.g.

Plot[CoshIntegral[x] Sinh[x] - Cosh[x] SinhIntegral[x], {x, 0, 40}, 
 WorkingPrecision -> 40, PlotRange -> All]

I verified this behavior. I tried to increase $MaxExtraPrecision but I didn't get any improvement.

A solution is to increase WorkingPrecision in Plot:

Plot[CoshIntegral[x] Sinh[x] - Cosh[x] SinhIntegral[x], {x, 0, 30}, 
 WorkingPrecision -> 50, PlotRange -> {{0, 30}, {-1, 1}}]

At last I created Cosh and Sinh integrals below with increased WorkingPrecision inside NIntegrate.

chi[z_] := N[EulerGamma + Log[z] + 
             NIntegrate[(Cosh[t] - 1)/t, {t, 0, z}, WorkingPrecision -> 50], 50]

shi[z_] := NIntegrate[Sinh[t]/t, {t, 0, z}, WorkingPrecision -> 50]

Then we can define res[x] as:

res[x_] := N[chi[x] Sinh[x] - shi[x] Cosh[x],50]

And Plot it:

Plot[res[x], {x, 0, 30}, 
 WorkingPrecision -> 50, PlotRange -> {{0, 30}, {-1, 1}}]]

It is a little bit slow but it works.

To check special values you can also use CoshIntegral Evaluation to verify the correctness for chi function above.

An alternative to setting the input precision is to ask for enough output accuracy. The main idea is to set the precision of the input to Infinity and use N to get the desired accuracy -- something like this:

N[f[SetPrecision[x, Infinity], 6]

Here's such an wrapper for a function:

Options[nEval] := {AccuracyGoal -> Automatic, "MaxExtraPrecision" -> Automatic};
nEval[f_, x_?NumericQ, OptionsPattern[]] := 
  Block[{$MaxExtraPrecision = 
     OptionValue["MaxExtraPrecision"] /. Automatic -> $MaxExtraPrecision},
   N[f[SetPrecision[x, Infinity]], OptionValue[AccuracyGoal] /. Automatic -> 6]
   ];

Then

myF[x_] := CoshIntegral[x] Sinh[x] - Cosh[x] SinhIntegral[x];

Plot[nEval[myF, x], {x, 0, 30}, PlotRange -> {{0, 30}, {-1, 1}}]

The precision is adapted to what is needed.

Plot[nEval[myF, x], {x, 30, 60}, PlotRange -> {Automatic, {-0.35, 0}}]

But one can still get a catastrophic loss of precision, due to the limit $MaxExtraPrecision.

Plot[nEval[myF, x], {x, 30, 100}, PlotRange -> {Automatic, {-0.35, 0}}]

In that case, one can raise the amount of extra precision allowed, even to Infinity. (In some cases, it might take an exorbitant amount of time to finish, but in this case, it is still quick.)

Plot[nEval[myF, x, "MaxExtraPrecision" -> Infinity], {x, 30, 100}, 
 PlotRange -> {Automatic, {-0.35, 0}}]