Some Strange results from MMA and Matlab

Recently I am repeating a numerical discrete Fourier transform.
Since this problem is related to a real experiment.

So the authors use some complex functions to fit the results.

Below formula is the Laplace transform of the waiting time PDF:

PHI30EQ[u_] := (u - 52/15*u^2 - 16/15*u^3)*Exp[u]*
    BesselK[1, u] + (7*u + 4*u^2 + 16/15*u^3)*Exp[u]*
    BesselK[0, u] - (2*u/Pi)^(1/2)*16/5;

Table[N[PHI30EQ[u], 30], {u, 0, 0.1, 0.001}]

The aim is to calculate the above formula numerically. The following information can be ignored.

The results given by MMA are:

Indeterminate,0.966028,0.96976,0.977704,0.987383,0.997947,1.00901,1.02035,1.03185,1.04345,1.05507,1.0667,1.0783,1.08987,1.10138,1.11284,1.12423,1.13555,1.1468,1.15798,1.16908,1.18011,1.19106,1.20193,1.21274,1.22346,1.23412,1.2447,1.25521,1.26565,1.27602,1.28632,1.29656,1.30672,1.31683,1.32687,1.33684,1.34676,1.35661,1.3664,1.37614,1.38581,1.39543,1.405,1.41451,1.42396,1.43336,1.44271,1.45201,1.46126,1.47046,1.47961,1.48871,1.49776,1.50677,1.51573,1.52464,1.53351,1.54234,1.55112,1.55986,1.56856,1.57721,1.58583,1.59441,1.60294,1.61144,1.6199,1.62832,1.6367,1.64504,1.65335,1.66162,1.66986,1.67806,1.68623,1.69436,1.70246,1.71052,1.71856,1.72656,1.73452,1.74246,1.75036,1.75824,1.76608,1.77389,1.78167,1.78942,1.79715,1.80484,1.81251,1.82014,1.82775,1.83533,1.84289,1.85041,1.85791,1.86539,1.87283,1.88026

Do you think the above numerical result is correct?

Table[N[PHI30EQ[u], 30], {u, 0, Rationalize[0.1], Rationalize[0.001]}]

Limit[PHI30EQ[u], u -> 0]

PHI30EQ[_?(EqualTo[0])] = 1;