f = 1.3 10^9; (*Hz*)
[Omega] = 2. Pi f;
Trf = 1./f;

Tc = 9.25; (*K*)
T = 1.5; (*K*)

[Mu]0 = 4. Pi 10^-7; (*H/m*)
[Phi]0 = 2.07 10^-15; (*Wb*)
n = 5.56 10^28; (*m^-3*)
e = 1.6 10^-19; (*C*)
m = 9.1 10^-31; (*kg*)

[Tau] = 0.5 10^-9; (*s*)
vf = 1.37 10^6; (*m/s*)
Bc2 = 410. 10^-3; (*T*)
l = 100 10^-9; (*m*)
[Lambda] = 49.5 10^-9; (*m*)
[Xi] = 28 10^-9; (*m*)
[Kappa] = [Lambda]/[Xi];

[Gamma] = ([Phi]0)/([Mu]0 [Lambda] );

v0 = Sqrt[(l vf Sqrt[14 Zeta[3] (1 - T/Tc)])/(3 Pi [Tau])];

[Rho] = (m vf)/(n e^2 l);
[Eta] = ([Phi]0 Bc2)/[Rho];

g = Log[[Kappa]]+0.5+Exp[-0.4-0.8 Log[[Kappa]]-0.1(Log[[Kappa]])^2];
[Epsilon] = ([Phi]0^2 g)/(4. Pi [Mu]0 [Lambda]^2);

Tmax = 3. Trf;
Zmax = 1. 10^-6; (*m*)

B = 100. 10^-3; (*T*)

sol = NDSolve[{
     (v0^2 [Eta])/(v0^2 + D[x[t, z], t]^2)D[x[t, z],t] == [Epsilon] D[x[t, z], z, z] + [Gamma] B Cos[[Omega] t] Exp[-(z/[Lambda])],
     x[0., z] == 0.,
     x[t, Zmax] == 0.,
     (D[x[t, z], z] /. z -> 0) == 0.
     }, x, {t, 0., Tmax}, {z, 0., Zmax},
    AccuracyGoal -> 13,
    PrecisionGoal -> 2,
    MaxSteps -> Infinity];
u[t_, z_] = Evaluate[x[t, z] /. sol][[1]];
du[t_, z_] = D[u[t, z],t];

phaseSpace =
   ParallelTable[{Re[u[t, 0.]] 10^6, Re[du[t, 0.]] 10^-3}, {t, Trf,
     Tmax, Tmax/500.}];
ListPlot[phaseSpace, Joined -> True, PlotRange -> All]