Needs["CCompilerDriver`"];
Needs["SymbolicC`"];

Data

Input parameters

SetDirectory[NotebookDirectory[]];

"Length Element";
ll=1;
"Number of discretization";
e=1;
"define factor m for wavelength λ";
m=1;
λ=ll/2 m;
"Material parameters";
μ=1; (*Shear moduli*)
NN =0.4;(*Hardening exponent*)
kk=0.87;(*Prestress*)
ec = FindRoot[NN == ϵc Tanh[ϵc], {ϵc, 0.1}];(*Prestrain*)
ϵc = ϵc /. ec;
ϵ[δ_] := (δ) ϵc;
ee = ArcTanh[kk]/2;
ee/ϵc;
u = NN/(2 ee Coth[2 ee]);(*ξ in the article, ration μ*/μ *)
θ=Pi/2 - ArcTan[Exp[ϵ[1]]];(*Shear band inclination*)
αw=θ;(*direction of wave*)
t=0.59; (*Factor for evaluation of η*)

Green’s functions/tractions

Ω=2* Pi/λ Sqrt[( μ(1+kk) )/ρ] ;(*Frequency*)
η=t*kk;
μs=u μ; (*Shear moduli*)
ρ=1;
σ2=η μ-kk μ;
σ1=2 kk μ+σ2;
p=(σ1+σ2)/2;
σ=σ1-σ2;
kk=σ/(2 μ);
Pi/m;

p=(σ1+σ2)/2;     
σ=σ1-σ2;
kk=σ/(2 μ);
Δ=kk^2-4 u+4 u^2;
γ1=(1-2 u+Sqrt[Δ])/(1+kk);
γ2=(1-2 u-Sqrt[Δ])/(1+kk);
Λ[ϕ_]:=Sin[ϕ]^4 (Cot[ϕ]^2-γ1)(Cot[ϕ]^2-γ2);
w1[ϕ_]:=Cos[ϕ];
w2[ϕ_]:=Sin[ϕ];
L[ϕ_]:=μ(1+kk) Λ[ϕ];
wx[α_]:= -Cos[α]^2 ;
Ξ[r_,φ_]:=Cos[r φ]CosIntegral[Abs[r φ]]+Sin[r φ] SinIntegral[r φ]-I Pi/2 Cos[r φ];
φ[α_]:=(Ω  ρ )/Sqrt[L[α+θ]] Cos[α];



t111[ϕ_]:=((2 μs-p) w1[ϕ]w2[ϕ]^2 +w1[ϕ]((2 μs-μ)w2[ϕ]^2+(μ+σ/2) w1[ϕ]^2)) 1/L[ϕ];
t121[ϕ_]:=-((μ-p) w1[ϕ]^2 w2[ϕ]-(μ-σ/2) w2[ϕ]^3) 1/L[ϕ];
t222[ϕ_]:=((2 μs-p) w1[ϕ]^2 w2[ϕ] +w2[ϕ]((2 μs-μ)w1[ϕ]^2+(μ-σ/2) w2[ϕ]^2)) 1/L[ϕ];
t212[ϕ_]:=-((μ-p) w1[ϕ] w2[ϕ]^2-(μ+σ/2) w1[ϕ]^3) 1/L[ϕ];
t122[ϕ_]:=(-(2 μs-p) w1[ϕ]w2[ϕ]^2 +w1[ϕ]((2 μs-μ)w2[ϕ]^2+(μ+σ/2) w1[ϕ]^2)) 1/L[ϕ];
t112[ϕ_]:=((μ-p) w2[ϕ]^3 -(μ+σ/2) w1[ϕ]^2 w2[ϕ]) 1/L[ϕ];
t211[ϕ_]:=(-(2 μs-p) w1[ϕ]^2 w2[ϕ] +w2[ϕ]((2 μs-μ)w1[ϕ]^2+(μ-σ/2) w2[ϕ]^2)) 1/L[ϕ];
t221[ϕ_]:=((μ-p) w1[ϕ]^3 -(μ-σ/2) w2[ϕ]^2 w1[ϕ]) 1/L[ϕ];




tc1111[ϕ_]:=t111[ϕ] w1[ϕ];
tc1112[ϕ_]:=t111[ϕ] w2[ϕ];
tc1211[ϕ_]:=t121[ϕ] w1[ϕ];
tc1212[ϕ_]:=t121[ϕ] w2[ϕ];
tc2221[ϕ_]:=t222[ϕ] w1[ϕ];
tc2222[ϕ_]:=t222[ϕ] w2[ϕ];
tc2121[ϕ_]:=t212[ϕ] w1[ϕ];
tc2122[ϕ_]:=t212[ϕ] w2[ϕ];
tc1221[ϕ_]:=t122[ϕ] w1[ϕ];
tc1222[ϕ_]:=t122[ϕ] w2[ϕ];
tc1121[ϕ_]:=t112[ϕ] w1[ϕ];
tc1122[ϕ_]:=t112[ϕ] w2[ϕ];
tc2111[ϕ_]:=t211[ϕ] w1[ϕ];
tc2112[ϕ_]:=t211[ϕ] w2[ϕ];
tc2211[ϕ_]:=t221[ϕ] w1[ϕ];
tc2212[ϕ_]:=t221[ϕ] w2[ϕ];


t1111[α_]:=(tc1111[α+θ]-tc1111[θ+Pi/2]-tc1111'[θ+Pi/2](α-Pi/2)) (1/wx[α]);
t1112[α_]:=(tc1112[α+θ]-tc1112[θ+Pi/2]-tc1112'[θ+Pi/2](α-Pi/2)) (1/wx[α]);
t1211[α_]:=(tc1211[α+θ]-tc1211[θ+Pi/2]-tc1211'[θ+Pi/2](α-Pi/2)) (1/wx[α]);
t1212[α_]:=(tc1212[α+θ]-tc1212[θ+Pi/2]-tc1212'[θ+Pi/2](α-Pi/2)) (1/wx[α]);
t2221[α_]:=(tc2221[α+θ]-tc2221[θ+Pi/2]-tc2221'[θ+Pi/2](α-Pi/2)) (1/wx[α]);
t2222[α_]:=(tc2222[α+θ]-tc2222[θ+Pi/2]-tc2222'[θ+Pi/2](α-Pi/2)) (1/wx[α]);
t2121[α_]:=(tc2121[α+θ]-tc2121[θ+Pi/2]-tc2121'[θ+Pi/2](α-Pi/2)) (1/wx[α]);
t2122[α_]:=(tc2122[α+θ]-tc2122[θ+Pi/2]-tc2122'[θ+Pi/2](α-Pi/2)) (1/wx[α]);
t1221[α_]:=(tc1221[α+θ]-tc1221[θ+Pi/2]-tc1221'[θ+Pi/2](α-Pi/2)) (1/wx[α]);
t1222[α_]:=(tc1222[α+θ]-tc1222[θ+Pi/2]-tc1222'[θ+Pi/2](α-Pi/2)) (1/wx[α]);
t1121[α_]:=(tc1121[α+θ]-tc1121[θ+Pi/2]-tc1121'[θ+Pi/2](α-Pi/2)) (1/wx[α]);
t1122[α_]:=(tc1122[α+θ]-tc1122[θ+Pi/2]-tc1122'[θ+Pi/2](α-Pi/2)) (1/wx[α]);
t2111[α_]:=(tc2111[α+θ]-tc2111[θ+Pi/2]-tc2111'[θ+Pi/2](α-Pi/2)) (1/wx[α]);
t2112[α_]:=(tc2112[α+θ]-tc2112[θ+Pi/2]-tc2112'[θ+Pi/2](α-Pi/2)) (1/wx[α]);
t2211[α_]:=(tc2211[α+θ]-tc2211[θ+Pi/2]-tc2211'[θ+Pi/2](α-Pi/2)) (1/wx[α]);
t2212[α_]:=(tc2212[α+θ]-tc2212[θ+Pi/2]-tc2212'[θ+Pi/2](α-Pi/2)) (1/wx[α]);


Off[NIntegrate::slwcon];

Off[NIntegrate::ncvb];
Off[General::stop];
ts1111=-2/(4 Pi^2) NIntegrate[t1111[α],{α,0,Pi},MaxRecursion->4];
ts1112=-2/(4 Pi^2) NIntegrate[t1112[α],{α,0,Pi},MaxRecursion->4];
ts1211=-2/(4 Pi^2) NIntegrate[t1211[α],{α,0,Pi},MaxRecursion->4];
ts1212=-2/(4 Pi^2) NIntegrate[t1212[α],{α,0,Pi},MaxRecursion->4];
ts2221=-2/(4 Pi^2) NIntegrate[t2221[α],{α,0,Pi},MaxRecursion->4];
ts2222=-2/(4 Pi^2) NIntegrate[t2222[α],{α,0,Pi},MaxRecursion->4];
ts2121=-2/(4 Pi^2) NIntegrate[t2121[α],{α,0,Pi},MaxRecursion->4];
ts2122=-2/(4 Pi^2) NIntegrate[t2122[α],{α,0,Pi},MaxRecursion->4];
ts1221=-2/(4 Pi^2) NIntegrate[t1221[α],{α,0,Pi},MaxRecursion->4];
ts1222=-2/(4 Pi^2) NIntegrate[t1222[α],{α,0,Pi},MaxRecursion->4];
ts1121=-2/(4 Pi^2) NIntegrate[t1121[α],{α,0,Pi},MaxRecursion->4];
ts1122=-2/(4 Pi^2) NIntegrate[t1122[α],{α,0,Pi},MaxRecursion->4];
ts2111=-2/(4 Pi^2) NIntegrate[t2111[α],{α,0,Pi},MaxRecursion->4];
ts2112=-2/(4 Pi^2) NIntegrate[t2112[α],{α,0,Pi},MaxRecursion->4];
ts2211=-2/(4 Pi^2) NIntegrate[t2211[α],{α,0,Pi},MaxRecursion->4];
ts2212=-2/(4 Pi^2) NIntegrate[t2212[α],{α,0,Pi},MaxRecursion->4];

TTst={
{{{ts1111(* k=1 *),ts1112(* k=2 *)}(* j=1 *),{ts1121(* k=1 *),ts1122(* k=2 *)}(* j=2 *)}(* i=1 *),{{ts1211(* k=1 *),ts1212(* k=2 *)}(* j=1 *),{ts1221(* k=1 *),ts1222(* k=2 *)}(* j=2 *)}(* i=2 *)},  (* g=1 *)

{{{ts2111(* k=1 *),ts2112(* k=2 *)}(* j=1 *),{ts2121(* k=1 *),ts2122(* k=2 *)}(* j=2 *)}(* i=1 *),{{ts2211(* k=1 *),ts2212(* k=2 *)}(* j=1 *),{ts2221(* k=1 *),ts2222(* k=2 *)}(* j=2 *)}(* i=2 *)} (* g=2 *)
};
Chop[MatrixForm[TTst]];

td1111[r_,α_]:=-2/(4 Pi^2) (tc1111[α+θ]-w1[α+θ]^2)(-(Ω^2/L[α+θ])Ξ[r,φ[α]]);
td1112[r_,α_]:=-2/(4 Pi^2) (tc1112[α+θ]-w1[α+θ] w2[α+θ])(-(Ω^2/L[α+θ])Ξ[r,φ[α]]);
td1211[r_,α_]:=-2/(4 Pi^2) (tc1211[α+θ])(-(Ω^2/L[α+θ])Ξ[r,φ[α]]);
td1212[r_,α_]:=-2/(4 Pi^2) (tc1212[α+θ])(-(Ω^2/L[α+θ])Ξ[r,φ[α]]);
td2221[r_,α_]:=-2/(4 Pi^2) (tc2221[α+θ]-w2[α+θ]w1[α+θ])(-(Ω^2/L[α+θ])Ξ[r,φ[α]]);
td2222[r_,α_]:=-2/(4 Pi^2) (tc2222[α+θ]- w2[α+θ]^2)(-(Ω^2/L[α+θ])Ξ[r,φ[α]]);
td2121[r_,α_]:=-2/(4 Pi^2) (tc2121[α+θ])(-(Ω^2/L[α+θ])Ξ[r,φ[α]]);
td2122[r_,α_]:=-2/(4 Pi^2) (tc2122[α+θ])(-(Ω^2/L[α+θ])Ξ[r,φ[α]]);
td1221[r_,α_]:=-2/(4 Pi^2) (tc1221[α+θ]-w1[α+θ]^2)(-(Ω^2/L[α+θ])Ξ[r,φ[α]]);
td1222[r_,α_]:=-2/(4 Pi^2) (tc1222[α+θ]-w1[α+θ] w2[α+θ])(-(Ω^2/L[α+θ])Ξ[r,φ[α]]);
td1121[r_,α_]:=-2/(4 Pi^2) (tc1121[α+θ])(-(Ω^2/L[α+θ])Ξ[r,φ[α]]);
td1122[r_,α_]:=-2/(4 Pi^2) (tc1122[α+θ])(-(Ω^2/L[α+θ])Ξ[r,φ[α]]);
td2111[r_,α_]:=-2/(4 Pi^2) (tc2111[α+θ]-w1[α+θ] w2[α+θ])(-(Ω^2/L[α+θ])Ξ[r,φ[α]]);
td2112[r_,α_]:=-2/(4 Pi^2) (tc2112[α+θ]-w2[α+θ]^2)(-(Ω^2/L[α+θ])Ξ[r,φ[α]]);
td2211[r_,α_]:=-2/(4 Pi^2) (tc2211[α+θ])(-(Ω^2/L[α+θ])Ξ[r,φ[α]]);
td2212[r_,α_]:=-2/(4 Pi^2) (tc2212[α+θ])(-(Ω^2/L[α+θ])Ξ[r,φ[α]]);

TTd[r_,α_]:={
{{{td1111[r,α](* k=1 *),td1112[r,α](* k=2 *)}(* j=1 *),{td1121[r,α](* k=1 *),td1122[r,α](* k=2 *)}(* j=2 *)}(* i=1 *),{{td1211[r,α](* k=1 *),td1212[r,α](* k=2 *)}(* j=1 *),{td1221[r,α](* k=1 *),td1222[r,α](* k=2 *)}(* j=2 *)}(* i=2 *)},  (* g=1 *)

{{{td2111[r,α](* k=1 *),td2112[r,α](* k=2 *)}(* j=1 *),{td2121[r,α](* k=1 *),td2122[r,α](* k=2 *)}(* j=2 *)}(* i=1 *),{{td2211[r,α](* k=1 *),td2212[r,α](* k=2 *)}(* j=1 *),{td2221[r,α](* k=1 *),td2222[r,α](* k=2 *)}(* j=2 *)}(* i=2 *)} (* g=2 *)
};


K={
{{{K1111(* k=1 *),K1112(* k=2 *)}(* h=1 *),{K1121(* k=1 *),K1122(* k=2 *)}(* h=2 *)}(* j=1 *),{{K1211(* k=1 *),K1212(* k=2 *)}(* h=1 *),{K1221(* k=1 *),K1222(* k=2 *)}(* h=2 *)}(* j=2 *)},  (* i=1 *)
{{{K2111(* k=1 *),K2112(* k=2 *)}(* h=1 *),{K2121(* k=1 *),K2122(* k=2 *)}(* h=2 *)}(* j=1 *),{{K2211(* k=1 *),K2212(* k=2 *)}(* h=1 *),{K2221(* k=1 *),K2222(* k=2 *)}(* h=2 *)}(* j=2 *)} (* i=2 *)
};
K1111:=μs-σ/2-p;
K1122:=-μs;
K2211:=-μs;
K2222:=μs+σ/2-p;
K1212:=μ+σ/2;
K2121:=μ-σ/2;
K1221:=μ-p;
K2112:=K1221;
K1112:=0;
K1121:=0;
K1211:=0;
K2111:=0;
K2221:=0;
K2212:=0;
K2122:=0;
K1222:=0;


n={-Sin[θ],Cos[θ]};
s={Cos[θ],Sin[θ]};
Q={{Cos[θ],-Sin[θ]},{Sin[θ],Cos[θ]}};

Tst[x_] = Chop[Sum[Sum[Sum[Sum[Sum[Sum[n[[l]] s[[mm]] K[[l,mm,k,g]] TTst[[g,ii,jj,k]] n[[ii]] s[[jj]] 1/Abs[x]^2,{l,1,2}],{mm,1,2}],{k,1,2}],{g,1,2}],{ii,1,2}],{jj,1,2}]];
cα=Tst[x] Abs[x]^2;
Td[r_, α_] = Chop[Sum[Sum[Sum[Sum[Sum[Sum[n[[l]] s[[mm]] K[[l,mm,k,g]] TTd[r,α][[g,ii,jj,k]] n[[ii]] s[[jj]],{l,1,2}],{mm,1,2}],{k,1,2}],{g,1,2}],{ii,1,2}],{jj,1,2}]];

Shape functions

QN1sx[Tα_]:=Tα ((1/(2le))-(3/(2le))Log[2 le]);
QN1c[Tα_]:=Tα (1/le);
QN1dx[Tα_]:=Tα ((1/le)-(1/(2le))Log[2 le]);

QN2sx[Tα_]:=Tα (-(2/le)+(2/le)Log[2 le]);
QN2c[Tα_]:=-Tα (4/le);
QN2dx[Tα_]:=Tα (-(2/le)+(2/le)Log[2 le]);

QN3sx[Tα_]:=Tα ((1/le)-(1/(2le))Log[2 le]);
QN3c[Tα_]:=Tα (1/le);
QN3dx[Tα_]:=Tα ((1/(2le))-(3/(2le))Log[2 le]);
ξN1QQ[r_,α_,j_]:=(r/(2le))-(((j-α))/2);
ξN2QQ[r_,α_,j_]:=(r/(2le))-(((j-α-1))/2);
ξN3QQ[r_,α_,j_]:=(r/(2le))-(((j-α-2))/2);


QN1GenSt[T_,ξ_,α_,j_]:=NIntegrate[T(1-3ξ+2(ξ^2)),{r,le (j-α),le (j+2-α)}];
QN2GenSt[T_,ξ_,α_,j_]:=NIntegrate[T( 4ξ-4(ξ^2)),{r,le (j-α-1),le (j+1-α)}];
QN3GenSt[T_,ξ_,α_,j_]:=NIntegrate[T(-ξ+2(ξ^2)),{r,le (j-α-2),le (j-α)}];

QN1GenD[T_,ξ_,y_,j_]:=NIntegrate[T(1-3ξ+2(ξ^2)),{α,0,Pi},{r,le (j-y),le (j+2-y)},MaxRecursion->4,Method->{"GlobalAdaptive","SymbolicProcessing"->0}];
QN2GenD[T_,ξ_,y_,j_]:=NIntegrate[T( 4ξ-4(ξ^2)),{α,0,Pi},{r,le (j-y-1),le (j+1-y)},MaxRecursion->4,Method->{"GlobalAdaptive","SymbolicProcessing"->0}];
QN3GenD[T_,ξ_,y_,j_]:=NIntegrate[T(-ξ+2(ξ^2)),{α,0,Pi},{r,le (j-y-2),le (j-y)},MaxRecursion->4,Method->{"GlobalAdaptive","SymbolicProcessing"->0}];

QPN1sx[Tα_]:=Tα (-(7/(16le))-(3/(4le))Log[2 le]);
QPN1c[Tα_]:=Tα ((1/(2le))+(1/(2le))Log[3]);
QPN1dx[Tα_]:=0;

QPN2sx[Tα_]:=0;
QPN2c[Tα_]:=-Tα (14/(3le));
QPN2dx[Tα_]:=Tα (-(1/(4le))+(1/le)Log[2 le]);

QPN3sx[Tα_]:=0;
QPN3c[Tα_]:=Tα ((1/(2le))+(1/(2le))Log[3]);
QPN3dx[Tα_]:=Tα (-(7/(16le))-(3/(4le))Log[2 le]);


QPN2GenStsx[T_,ξ_,α_,j_]:=NIntegrate[T( 4(Sqrt[ξ])-4ξ),{r,le (j-α-1),le (j+1-α)}];
QPN3GenStsx[T_,ξ_,α_,j_]:=NIntegrate[T(-(Sqrt[ξ])+2ξ),{r,le (j-α-2),le (j-α)}];

QPN1GenStdx[T_,ξ_,α_,j_]:=NIntegrate[T(-(Sqrt[1-ξ])-2(-1+ξ)),{r,le (j-α),le (j-α+2)}];
QPN2GenStdx[T_,ξ_,α_,j_]:=NIntegrate[T( 4(Sqrt[1-ξ])+4(-1+ξ)),{r,le (j-α-1),le (j+1-α)}];


QPN2GenDsx[T_,ξ_,y_,j_]:=NIntegrate[T( 4(Sqrt[ξ])-4ξ),{α,0,Pi},{r,le (j-y-1),le (j+1-y)},MaxRecursion->4,Method->{"GlobalAdaptive","SymbolicProcessing"->0}];
QPN3GenDsx[T_,ξ_,y_,j_]:=NIntegrate[T(-(Sqrt[ξ])+2ξ),{α,0,Pi},{r,le (j-y-2),le (j-y)},MaxRecursion->4,Method->{"GlobalAdaptive","SymbolicProcessing"->0}];

QPN1GenDdx[T_,ξ_,y_,j_]:=NIntegrate[T(-(Sqrt[1-ξ])-2(-1+ξ)),{α,0,Pi},{r,le (j-y),le (j-y+2)},MaxRecursion->4,Method->{"GlobalAdaptive","SymbolicProcessing"->0}];
QPN2GenDdx[T_,ξ_,y_,j_]:=NIntegrate[T( 4(Sqrt[1-ξ])+4(-1+ξ)),{α,0,Pi},{r,le (j-y-1),le (j+1-y)},MaxRecursion->4,Method->{"GlobalAdaptive","SymbolicProcessing"->0}];

Statics

Ast:=Table[a[i,j],{i,e-1},{j,e-1}];
le=ll/e;
A11=ParallelTable[a1[α,α]=QN2c[cα],{α,3,e-3,2}];AppendTo[A11,QPN2c[cα]];PrependTo[A11,QPN2c[cα]];
A22=ParallelTable[a1[α,α]=QN1sx[cα]+QN3dx[cα],{α,4,e-4,2}];

AppendTo[A22,QN1sx[cα]+QPN3dx[cα]];
PrependTo[A22,QN3dx[cα]+QPN1sx[cα]];

A12=ParallelTable[j=α+1;
a1[α,α+1]=QN3c[cα]+QN1GenSt[Tst[r],ξN1QQ[r,α,j],α,j],{α,3,e-5,2}];Clear@j;
PrependTo[A12,QPN3c[cα]+QN1GenSt[Tst[r],ξN1QQ[r,0.5,2],0.5,2]];
AppendTo[A12,QN3c[cα]+QPN1GenStdx[Tst[r],ξN1QQ[r,e-3,e-2],e-3,e-2]];

A23=ParallelTable[a1[α,α+1]=QN2sx[cα],{α,2,e-4,2}];AppendTo[A23,QPN2dx[cα]];

A24=ParallelTable[j=α+2;a1[α,α+2]=QN3sx[cα]+QN1GenSt[Tst[r],ξN1QQ[r,α,j],α,j],{α,2,e-6,2}];Clear@j;
AppendTo[A24,α=e-4;j=α+2;QN3sx[cα]+QPN1GenStdx[Tst[r],ξN1QQ[r,α,j],α,j]];
Adisp1=ParallelTable[a1[α,j]=QN2GenSt[Tst[r],ξN2QQ[r,α,j],α,j],{α,3,e-5,2},{j,α+2,e-3,2}];
PrependTo[Adisp1,ParallelTable[a1[1,j]=QN2GenSt[Tst[r],ξN2QQ[r,0.5,j],0.5,j],{j,3,e-3,2}]];
Clear@α
Adisp2=ParallelTable[a1[α,j]=QN2GenSt[Tst[r],ξN2QQ[r,α,j],α,j],{α,2,e-6,2},{j,α+3,e-3,2}];

Adisp3=ParallelTable[a1[α,j]=QPN2GenStdx[Tst[r],ξN2QQ[r,α,j],α,j],{α,2,e-3,1},{j,e-1,e-1}];
PrependTo[Adisp3,ParallelTable[QPN2GenStdx[Tst[r],ξN2QQ[r,0.5,j],0.5,j],{j,e-1,e-1}]];
Clear@α


Apari1=ParallelTable[a1[α,j]=QN1GenSt[Tst[r],ξN1QQ[r,α,j],α,j]+QN3GenSt[Tst[r],ξN3QQ[r,α,j],α,j],{α,3,e-7,2},{j,α+3,e-4,2}];
PrependTo[Apari1,ParallelTable[a1[1,j]=QN1GenSt[Tst[r],ξN1QQ[r,0.5,j],0.5,j]+QN3GenSt[Tst[r],ξN3QQ[r,0.5,j],0.5,j],{j,4,e-4,2}]];
Clear@α
Apari2=ParallelTable[a1[α,j]=QN1GenSt[Tst[r],ξN1QQ[r,α,j],α,j]+QN3GenSt[Tst[r],ξN3QQ[r,α,j],α,j],{α,2,e-8,2},{j,α+4,e-4,2}];

Apari3=ParallelTable[a1[α,j]=QN3GenSt[Tst[r],ξN3QQ[r,α,j],α,j]+QPN1GenStdx[Tst[r],ξN1QQ[r,α,j],α,j],{α,2,e-5,1},{j,e-2,e-2}];
PrependTo[Apari3,ParallelTable[QN3GenSt[Tst[r],ξN3QQ[r,0.5,j],0.5,j]+QPN1GenStdx[Tst[r],ξN1QQ[r,0.5,j],0.5,j],{j,e-2,e-2}]];
Clear@α

For[i=1,i<(e-1)/2+1,i++,a[2i-1,2i-1]=A11[[i]]];
For[i=1,i<(e-1)/2,i++,a[2i,2i]=A22[[i]]];
For[i=1,i<(e-1)/2,i++,a[2i-1,2i]=A12[[i]]];
For[i=1,i<(e-1)/2,i++,a[2i,2i+1]=A23[[i]]];
For[i=1,i<(e)/2-1,i++,a[2i,2i+2]=A24[[i]]];


For[i=1,i<(e)/2-1,i++,For[j=1,j<(e)/2-i,j++,a[2i-1,2j+2i-1]=Adisp1[[i,j]]]];
For[i=1,i<(e)/2-2,i++,For[j=1,j<(e)/2-i-1,j++,a[2i,2j+2i+1]=Adisp2[[i,j]]]];
For[i=1,i<(e)/2-2,i++,For[j=1,j<(e)/2-i-1,j++,a[2i-1,2j+2i]=Apari1[[i,j]]]];
For[i=1,i<(e)/2-3,i++,For[j=1,j<(e)/2-i-2,j++,a[2i,2j+2i+2]=Apari2[[i,j]]]];
For[i=1,i<e-2,i++,a[i,e-1]=Adisp3[[i]][[1]]];

For[i=1,i<e-4,i++,a[i,e-2]=Apari3[[i]][[1]]];
For[i=1,i<e-1,i++,For[j=1,j<e-i,j++,a[e-i,e-i-j]=Ast[[i,j+i]]]];

Dynamics

Clear[α,j];

Ad:=Table[ad[i,j],{i,e-1},{j,e-1}];
le=ll/e;
Off[NIntegrate::slwcon];
Off[NIntegrate::nlim];
Off[NIntegrate::ncvb];
Off[General::stop];
D1=ParallelTable[Off[NIntegrate::slwcon];Off[NIntegrate::ncvb];Off[General::stop];d1[y,j]=QN2GenD[Td[r,α],ξN2QQ[r,y,j],y,j],{y,2,e-2},{j,3,e-3,2}];

PrependTo[D1,ParallelTable[Off[NIntegrate::slwcon];Off[NIntegrate::ncvb];Off[General::stop];y=0.5;d1[1,j]=QN2GenD[Td[r,α],ξN2QQ[r,y,j],y,j],{j,3,e-3,2}]]; Clear@y;
AppendTo[D1,ParallelTable[Off[NIntegrate::slwcon];Off[NIntegrate::ncvb];Off[General::stop];y=e-0.5;d1[e-1,j]=QN2GenD[Td[r,α],ξN2QQ[r,y,j],y,j],{j,3,e-3,2}]]; Clear@y;
D2=ParallelTable[Off[NIntegrate::slwcon];Off[NIntegrate::ncvb];Off[General::stop];d2[y,j]=QN1GenD[Td[r,α],ξN1QQ[r,y,j],y,j]+QN3GenD[Td[r,α],ξN3QQ[r,y,j],y,j],{y,2,e-2},{j,4,e-4,2}];
PrependTo[D2,ParallelTable[Off[NIntegrate::slwcon];Off[NIntegrate::ncvb];Off[General::stop];y=0.5;d2[1,j]=QN1GenD[Td[r,α],ξN1QQ[r,y,j],y,j]+QN3GenD[Td[r,α],ξN3QQ[r,y,j],y,j],{j,4,e-4,2}]]; Clear@y;
AppendTo[D2,ParallelTable[Off[NIntegrate::slwcon];Off[NIntegrate::ncvb];Off[General::stop];y=e-0.5;
d2[e-1,j]=QN1GenD[Td[r,α],ξN1QQ[r,y,j],y,j]+QN3GenD[Td[r,α],ξN3QQ[r,y,j],y,j],{j,4,e-4,2}]]; Clear@y;



j=1;D3=ParallelTable[Off[NIntegrate::slwcon];Off[NIntegrate::ncvb];Off[General::stop];d1[y,1]=QPN2GenDsx[Td[r,α],ξN2QQ[r,y,j],y,j],{y,2,e-2}];Clear@j;
PrependTo[D3,j=1;y=0.5;d1[1,1]=QPN2GenDsx[Td[r,α],ξN2QQ[r,y,j],y,j]];Clear[y,j];
AppendTo[D3,j=1;y=e-0.5;d1[e-1,1]=QPN2GenDsx[Td[r,α],ξN2QQ[r,y,j],y,j]];Clear[y,j];
D4=Reverse[D3];



D5=ParallelTable[Off[NIntegrate::slwcon];Off[NIntegrate::ncvb];Off[General::stop];
j=2;d2[y,2]=QPN3GenDsx[Td[r,α],ξN3QQ[r,y,j],y,j]+QN1GenD[Td[r,α],ξN1QQ[r,y,j],y,j],{y,2,e-2}];Clear[y,j];
PrependTo[D5,y=0.5;j=2;d2[1,2]=QPN3GenDsx[Td[r,α],ξN3QQ[r,y,j],y,j]+QN1GenD[Td[r,α],ξN1QQ[r,y,j],y,j]];Clear[y,j];
AppendTo[D5,y=e-0.5;j=2;d2[e-1,2]=QPN3GenDsx[Td[r,α],ξN3QQ[r,y,j],y,j]+QN1GenD[Td[r,α],ξN1QQ[r,y,j],y,j]];

D6=Reverse[D5];

For[i=1,i<e,i++,For[j=1,j<e/2-1,j++,ad[i,2j+1]=D1[[i,j]]]];
For[i=1,i<e,i++,For[j=1,j<e/2-2,j++,ad[i,2j+2]=D2[[i,j]]]];
For[i=1,i<e,i++,ad[i,1]=D3[[i]]];
For[i=1,i<e,i++,ad[i,e-1]=D4[[i]]];
For[i=1,i<e,i++,ad[i,2]=D5[[i]]];
For[i=1,i<e,i++,ad[i,e-2]=D6[[i]]];

Results

A=Ast+Ad;
Clear[i,F,V,f,v];
V:=Array[v,{e-1}];
F:=Array[f,{e-1}];
t={{0,0},{1,0}};
tt=Q.t.Q\[Transpose];
ttt=n.tt.s;


nn=αw;
λt=Sqrt[Λ[nn]] λ;
kw=2 Pi/λt;
τ=-1;
n1=Cos[nn];
n2=Sin[nn];

For[i=1,i<e,i++, f[i]=-τ/μ Exp[I kw (n1 (Cos[θ] (-ll/2+i le))+n2 (Sin[θ] (-ll/2+i le)))] ((n1^2 (p-μ)-(-1+kk) n2^2 μ) Cos[θ]^2+(n2^2 (p-μ)+(1+kk) n1^2 μ) Sin[θ]^2+n1 n2 (p-2 μs) Sin[2 θ])] ;
MatrixForm[A.V];
sol:=Solve[A.V==F,V];
sol;
sol=Flatten[sol];
X:=Array[x,{e-1}];

res2={{0,0}};res2=AppendTo[res2,{le/2,v[1]/.sol}];
For[i=2,i<e-1,i++,x[i]=i ll/e];
For[i=2,i<e-1,i++,res2=AppendTo[res2,{x[i],v[i]/.sol}]];
res2=AppendTo[res2,{ll-le/2,v[e-1]/.sol}];res2=AppendTo[res2,{ll,0}];
res2;

ListLinePlot[{Re@res2},Mesh->All]