At a glance one can see that the highest order derivatives for \[Mu]1 and P1 occur only in one equation, eqn6. Further, eqn5 involves only lower order derivatives (orders 1, 0, resp) of \[Mu]1, P1. So there are only two equations involving the three highest-order derivatives, eqn4 and eqn6. Consequently NDSolve "cannot solve to find an explicit formula for the [highest-order] derivatives." A simple fix is to differentiate eqn5. We also need to remove the initial value for P1'[0], since it is determined by the equations.
a1 = 1.38*10^10;
q1 = 1.6*10^-19;
\[Epsilon]1 = 8.85*10^-12;
\[Lambda]1 = 4*10^-9;
\[Mu]hat1 = 0.8*10^-20;
f1 = 0;
\[Beta]1 = 4.63*10^-44;
\[Alpha]1 = -80*10^-16;
M1 = 70*10^6;
eqn4 = (M1 - ((\[Alpha]1^2)/\[Beta]1)) u1''[
y] + (\[Alpha]1/\[Beta]1)*\[Mu]1'[y] == 0;
eqn5 = D[a1*P1[y] - (1/q1)*\[Mu]1'[y] == 0, y];
eqn6 = \[Mu]1''[y] + (q1/\[Epsilon]1)*
P1'[y] - (1/\[Lambda]1^2)*\[Mu]1[
y] + (1/\[Lambda]1^2)*\[Mu]hat1 == 0;
sol = NDSolve[{eqn4, eqn5,
eqn6, \[Mu]1'[0] == 0, \[Mu]1[0] == -8*10^-21, P1[0] == 10,(*P1'[
0]==0,*)u1[0] == 10^-9, u1'[0] == 0}, {\[Mu]1, P1, u1}, {y, 0,
10^-9}]
\[Mu]sol1[y_] := \[Mu]1[y] /. sol[[1]]
Psol1[y_] := P1[y] /. sol[[1]]
usol1[y_] := u1[y] /. sol[[1]]
