Stokes-Line Drawer for Landau-Zener type differential equations

1.Find Q[t] for the Landau-Zener type differential equation

Our starting point is
ι h (d/dt) (({{ψ1}, {ψ2}}))=(({{H11, H12}, {H21, H22}}))(({{ψ1}, {ψ2}})).
Here we assume
(({{H11, H12}, {H21, H22}}))=(({{k, (e^(i θ[t])) m0}, {(e^(-i θ[t])) m0, -k}}))
and decouple the equations to find Q in the standard form of the Exact WKB:
(((d^2) (h^2) ψ)/(dt^2))=Q(t) ψ

ClearAll["Global`*"];
(*===Assumptions===*)
H11=k;
H22=-k;
H12= m0 Exp[ I  θ[t]];
H21=  m0 Exp[ -I  θ[t]];
(*Show the Hamiltonian *)
"Your Hamiltonian is H=" MatrixForm[{{H11,H12},{H21,H22}}]
sol2=Part[Solve[I D[u1[t],t]==H11*u1[t]+H12*u2[t], u2[t]],1,1];
sol1=Part[Solve[I D[u2[t],t]==H21*u1[t]+H22*u2[t], u1[t]],1,1];
eq1=Simplify[I D[u2[t]/.sol2,t]-H21* u1[t]-H22*u2[t]/.sol2];
eq2=Simplify[I D[u1[t]/.sol1,t]-H12*u2[t]-H11*u1[t]/.sol1];
eq1normal = Simplify[eq1/Coefficient[eq1,u1''[t]]];
(* Check the coefficients of the decoupled equation u1'', u1', u1 for [c1 u1''+ c2 u1' + c3 u1=0] *)
"Decoupled equation for ψ1-> u1 is"  eq1normal
Coefficient[eq1normal,u1'[t]];
Coefficient[eq1normal,u1[t]];
eq2normal = Simplify[eq2/Coefficient[eq2,u2''[t]]];
"Decoupled equation for ψ2-> u2 is"  eq2normal
Coefficient[eq2normal,u2''[t]];
Coefficient[eq2normal,u2'[t]];
Coefficient[eq2normal,u2[t]];
(* Transform into the standard EWKB form*)
kata[t_]:=Simplify[Integrate[Coefficient[eq1normal,u1'[t]]/2,t]];
Simplify[D[Exp[-kata[t]]*hu1[t],t]];
heq10=eq1normal/.u1[t]->Exp[-kata[t]]*hu1[t];
heq11=heq10/.u1'[t]->D[Exp[-kata[t]]*hu1[t],t];
heq12=heq11/.u1''[t]->D[D[Exp[-kata[t]]*hu1[t],t],t];
heq10;
heq11;
heq12;
heq12normal=heq12/Coefficient[heq12,hu1''[t]];

Coefficient[heq12normal,hu1''[t]];
Simplify[Coefficient[heq12normal,hu1'[t]]];
Q1=Simplify[Coefficient[heq12normal,hu1[t]]];
"Q for ψ1(hat) is" Q1
kata2[t_]:=Simplify[Integrate[Coefficient[eq2normal,u2'[t]]/2,t]];
Simplify[D[Exp[-kata2[t]]*hu2[t],t]];
heq20=eq2normal/.u2[t]->Exp[-kata2[t]]*hu2[t];
heq21=heq20/.u2'[t]->D[Exp[-kata2[t]]*hu2[t],t];
heq22=heq21/.u2''[t]->D[D[Exp[-kata2[t]]*hu2[t],t],t];
heq20;
heq21;
heq22;
heq22normal=heq22/Coefficient[heq22,hu2''[t]];

Coefficient[heq12normal,hu2''[t]];
Simplify[Coefficient[heq12normal,hu2'[t]]];
Q2=Simplify[Coefficient[heq22normal,hu2[t]]];
"Q for ψ2(hat) is" Q2

θ[t_]:=A/(1+Exp[-a t])
a=1.;
A=Pi/3.;
m0=1;
k=1.0;
"For θ[t_]:=A/(1+Exp[-a t]) with[a=1.; A=Pi/3.; m0=1; k=1.0];"
"Your Q1 is"(Q1/.θ'[t]->D[θ[t],t])/.θ''[t]->D[D[θ[t],t],t]
"Your Q2 is"(Q2/.θ'[t]->D[θ[t],t])/.θ''[t]->D[D[θ[t],t],t]

2. Draw the stokes lines

Please put your Q manually in the code.

ClearAll["Global`*"];
(* FindTP(find turning points) and FindSP(find singular points) has to be customized for your Q *)
FindTP[Q_]:=Module[{x},
Flatten[x/.NSolve[Q[x]==0]/.(C[1])->{-1,0,1}]];
FindSP[Q_]:=Module[{x},
Flatten[x/.NSolve[1/Q[x]==0]/.(C[1])->{-1,0,1}]];
(* This function finds the initial three branches of the stokes lines at each Turning points *)
SolTPbranch[a_]:=Module[{dx,x,y,bp},
dx:=0.001(x+ I y);
bp=(x+ I y)/.NSolve[Im[Sqrt[Q'[a]dx]dx]==0&&{x,y}∈Circle[],{x,y}];
Arg[bp]];
(* StokesGetRep and Stokesfind1 finds the most probable direction from the given point on the stokes lines*)
StokesGetRep[xo_,θo_]:=Module[{ang,i},
ang[1]=θo;
For[i=2,i<iMax,i++,
ang[i]=Stokesfind1[xo,ang[i-1], Δθ PMi^(-i+2)]
];
ang[iMax-1]
];
Stokesfind1[xo_,θo_,dθ_]:=Module[{f,θ,z,mini,i,tab},
f[θ_]:=NIntegrate[Sqrt[Q[z]],{z,xo,N[xo +  dr (Cos[θ]+ I Sin[θ])]},WorkingPrecision->6];
 mini=Ordering[Table[Abs[Im[f[θo-dθ+ dθ i/PMi]]],{i,0,2 PMi}],1][[1]];
θo-dθ+ dθ (mini-1)/PMi];
(* GetStokesTable generates tables of the points on the stokes lines using StokesGetRep*)
GetStokesTable[k_,j_]:=Module[{i,nn},
sX[[k,j,1]]=TP[[k]];
angS[[k,j,1]]=TPbranch[[k,j]];
For[i=2,i<stepN+1,i++,
sX[[k,j,i]]=sX[[k,j,i-1]]+N[dr E^(I angS[[k,j,i-1]])];
SPTP=Join[TP,SP];
domain=Table[Re[SPTP[[nn]]]-0.04<Re[sX[[k,j,i]]]<Re[SPTP[[nn]]]+0.04&&Im[SPTP[[nn]]]-0.04<Im[sX[[k,j,i]]]<Im[SPTP[[nn]]]+0.04,{nn,nTP+nSP}];
If[i>10&&Length[DeleteDuplicates[domain]]==2,Break[]];
If[Abs[sX[[k,j,i]]]>Swindow,Break[]];(*縦幅に変更した*)
angS[[k,j,i]]= StokesGetRep[sX[[k,j,i]],angS[[k,j,i-1]]]
];
Table[{sX[[k,j,n]]},{n,1,stepN}]]

(* Your setups here*)
Q[t_]:= -(2.` +(0.2741556778080377` (E^(-2.` t)))/(((1+(E^(-1.` t))))^4)+(1.0471975511965976` (E^(-1.` t)))/(((1+(E^(-1.` t))))^2)+(1/2) I ((2.0943951023931953` (E^(-2.` t)))/(((1+(E^(-1.` t))))^3)-(1.0471975511965976` (E^(-1.` t)))/(((1+(E^(-1.` t))))^2)))
dr=0.01;(* Step of the move *)
stepN=500;(* Maximum number of steps for each stokes line *)
Δθ = 0.4;(* StokesGetRep searches direction in the range -Δθ-> +Δθ. Make this parameter small if stokes lines merge.*)
PMi=10;(* The search step of the direction is Δθ/PMi, repeated iMax times. *)
iMax=5;(* Depth of the direction search *)
Swindow=5.;(* StokesGetRep stops when Abs[sx]>Swindow or reaches at another turning points or singular points*) 
(*========本体===========================*)
TP=FindTP[Q];
SP=FindSP[Q];
"Your turning points are at" MatrixForm[TP]
"Your singular points are at" MatrixForm[SP]
nTP= Length[TP](* number of Turning Points *);
nSP= Length[SP](* number of Singular Points *);
TPbranch=Table[SolTPbranch[TP[[i]]],{i,nTP}];
MatrixForm[TPbranch];
(* Initializing the Matrices*)
bt=ConstantArray[N[0],{nTP,3}];
(*Table[0,{k,nTP},{j,3}];*)
sX=ConstantArray[0,{nTP,3,stepN}];
(*Table[0,{k,nTP},{j,3},{i,stepN}];*)
angS=ConstantArray[0,{nTP,3,stepN}];
(*Table[0,{k,nTP},{j,3},{i,stepN}];*)
AbsoluteTiming[StokesList=ParallelTable[Catenate[GetStokesTable[k,j]],{k,1,nTP},{j,1,3}];]
Show[ComplexListPlot[Catenate[StokesList],PlotRange->{-3-4I,3+4I}],ComplexListPlot[{TP,SP},PlotMarkers->{{●,50},{▲,50}},PlotRange->{-3-4I,3+4I}]]