I am trying to use NDSolve to find the numerical solution to a system of 4 first order PDEs with 8 boundary conditions (I think they are Dirichlet type conditions). There are 4 dependant variables and 2 independant variables. I am using cylindrical coordinates (Rr, Ttheta, Zz) but Ttheta is not used at all and so is excluded (may be the cuase of the problems?). I cannot for the life of me get Mathematica to run the code, I must be doing something wrong but I can’t seem to figure out what. If anyone is able to help in anyway, that would be much appreciated.

The error that this code creates changes quite often depending on little ‘fixes’ I have attempted to make, however I have been slowly learning that a lot of the error messages often don’t point to the cause of the problem so I have left them out of this post. The relevant equations in the code are also in the tradition form bellow the code just for ease of reading.

(*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] 

Thanks!