(*Constants*)
mi = 1.67*10^-27;
Vn = 10^4 ;
T = 5;
K = 8.617*10^-5;
Z = 1;

EQ0 = Vr[Rr, Zz]*mi*D[Log[ni[Rr, Zz]], Rr] +
   Vz[Rr, Zz]*mi*D[Log[ni[Rr, Zz]], Zz] + mi*D[Vr[Rr, Zz], Rr] +
   Vr[Rr, Zz]*mi + mi*D[Vz[Rr, Zz], Zz] == 0

EQ1 = -Vt[Rr,
     Zz]^2 == -Rr*((Z*K)/mi)*(D[T, Rr] + T*D[Log[ni[Rr, Zz]], Rr]) -
   Rr*Vn*Vr[Rr, Zz]

EQ2 = Vr[Rr, Zz]*D[Vt[Rr, Zz], Rr]*Rr + Vr[Rr, Zz]*Vt[Rr, Zz] == -Vn*
   Vt[Rr, Zz]*Rr

EQ3 = Vr[Rr, Zz]*D[Vz[Rr, Zz], Rr] +
   Vz[Rr, Zz]*
    D[Vz[Rr, Zz],
     Zz] == -((Z*K)/mi)*(D[T, Zz] + T*D[Log[ni[Rr, Zz]], Zz]) -
   Vn*Vz[Rr, Zz]


sol = NDSolve[{EQ0, EQ1, EQ2, EQ3,
    Vz[Rr, 0] == 10^3,
    Vt[Rr, 0] == 10^5*Rr,
    Vr[Rr, 0] == 0,
    Log[ni[Rr, 0]] == -(Rr^2/0.028^2) + Log[8.5*10^8],
    Vz[0.028, Zz] == 0,
    Vt[0, Zz] == 0,
    Vr[0, Zz] == 0,
    Log[ni[0.028, Zz]] == 0},
    {Vz, Vt, Vr, ni},
    {Rr, 0, 0.028}, {Zz, 0, 0.6}, AccuracyGoal -> 2,
  PrecisionGoal -> 2]