# Diffusion equation with concentration dependent diffusivity, NDSolve has problem at C0> 10^20

I am using NDSolve to solve the diffusion equation:

D[Cp[x, t], t] == D[Dc D[Cp[x, t], x], x]


where Dc is a function of Cp: Dc = Dc(Cp). The analytical expression for Dc is given. Here is my program.

Clear[DPi, DP, DIi, DI, DPI, DPVi, DPV, Dtotal, Tc, ni, n, Cp, nn, n1, n2, h]

q = 1.6 10^-19;kB = 1.38 10^-23;

ni[Tc_]:=1.8 10^21 Exp[-((0.66 q)/(kB Tc))];

n[Cp_, Tc_]:= 0.5 (Cp + Sqrt[Cp^2 + 4 ni[Tc]^2]);

h[Cp_, Tc_]:= 1 + n[Cp, Tc]/(2 ni[Tc]) ((n[Cp, Tc]/(2 ni[Tc]))^2 + 1)^-0.5;

DP[Cp_, Tc_]:=4.9 10^-5 Exp[-2.1 q/kB/Tc] (n[Cp, Tc]/ni[Tc])^0.5 h[Cp, Tc];

DI[Cp_, Tc_]:=0.23 Exp[-2.6 q/kB/Tc] (n[Cp, Tc]/ni[Tc])^-1.8 h[Cp, Tc];

DPI[Cp_, Tc_]:= (DP[Cp, Tc] DI[Cp, Tc])/(DP[Cp, Tc] + DI[Cp, Tc]);

DPV[Cp_, Tc_] :=7.6 10^4 Exp[-5.2 q/kB/Tc] (n[Cp, Tc]/ni[Tc])^2 h[Cp, Tc];

Dtotal[Cp_, Tc_] := DPI[Cp, Tc] + DPV[Cp, Tc];

Dc = Dtotal[Cp, 900 + 273] /. Cp -> Cp[x, t];

LogLogPlot[Dtotal[Cp, 850 + 273], {Cp, 10^15, 7 10^20},
PlotRange -> {{10^16, 10^21}, {10^-15, 10^-12}}, Frame -> True]

C0 = 3 10^20;

sol = NDSolve[{D[Cp[x, t], t] == D[ Dc D[Cp[x, t], x], x],
Cp[x, 0] == C0 UnitStep[-x], Cp[0, t] == C0,
Cp[3. 10^-4, t] == 0}, Cp, {x, 0, 3. 10^-4}, {t, 0, 5. 10^3},
Method -> {"MethodOfLines"}, MaxSteps -> 10000, PrecisionGoal -> 2];

LogPlot[Evaluate[Cp[x, 2400] /. sol], {x, 0, 1 10^-4},
PlotRange -> {{0, 1 10^-4}, {10^15, 1 10^21}}, Frame -> True ]


At low concentration from 10^18-10^20, I dont have any problem, but with concentrations higher than 10^20, the solution looks strange.