Mathematica does not evaluate the derivative of the following MeijerG function.
<< NumericalCalculus`;
ND[
MeijerG[{{2/3, 2/3, -(7/3), b1}, {13/6}}, {{5/3, 5/3, 7/6}, {5/3, 5/3}}, 1/5]
, {b1, 1}, 5/3]
Why is that?
Notice that the arguments are compatible with the usual rules for the MeijerG function, i.e.
ND[
MeijerG[{{2/3, 2/3, -(7/3), b1}, {13/6}}, {{5/3, 5/3, 7/6}, {5/3, 5/3}}, 1/5]
, {b1, 0}, 5/3]
makes sense and can be evaluated with Mathematica.
It seems that in evaluating the derivative ND is transforming the function to
MeijerG[{{-(7/3), 2/3, 2/3, 5/3 + 1}, {13/6}}, {{7/6, 5/3, 5/3}, {5/3, 5/3}}, 1/5]
which is ill defined. Perhaps there is such a rule that holds for symbolic derivatives but it's not defined for generic arguments. As it is right now, it looks like a bug and you could report it.
If you need just the numerical derivative, you can code it by hand via the definition. $$ \lim_{n\to\infty} n\big(f\big(x+\tfrac1n\big)-f(x)\big)\,. $$
If you want to get precise results you can even use LinearModelFit to estimate the error, like so
f[b1_] := MeijerG[{{2/3, 2/3, -(7/3), b1}, {13/6}}, {{5/3, 5/3, 7/6}, {5/3, 5/3}}, 1/5];
Table[{nn, (f[SetPrecision[5/3 + 1/nn, 50]] - f[SetPrecision[5/3, 50]]) nn}, {nn, 10^8, 10^9, 10^8}];
ListPlot[%]
fit = LinearModelFit[%%, {1}, x];
{fit["BestFit"], Max[Abs /@ fit["FitResiduals"]]}
(* {numerical derivative, error estimate} *)
This is probably like shooting a fly with a cannon, but it works.
Adjust the ranges of nn so that the plot looks approximately constant and adjust the precision inside SetPrecision so that the plot does not look too noisy.
b1have an assigned value? $\endgroup$