我的代码是个具有专业背景的问题，我将阐述背景情况：
此拉拔过程是复合金属丝的拉拔过程，就是有包层，金属丝内芯是材料B，包层材料是A。
因为A，B之间的流变应力阈值kf不一样，分为两种情况：A（Al99.5）小于B（SF-Cu），和A（SF-Cu）大于B（Al99.5)。
rA0 是整个复合金属丝横截面的半径3.5mm；
rB0 是金属丝中内芯横截面的半径3.1mm；
rA1 是拉拔完成后的整个金属丝的半径3mm；
拉拔口的拉角是 8度；
拉拔速度为0.01m/s；
还有A，B材料分别的比热容，密度和热传导能力已经给出
z 是拉拔方向上行进的距离，0是刚要进入形变区开始形变的位置。
开始的代码是数据调用这部分没有问题。还有DGLVec的公式集经老师查证都是对的，根据软件的错误显示错误应该是出在NDSolve命令中，我询问老师老师暂时顾不上我。无奈我已经1个月没有找到改进办法了！

以下是我的代码：
FileNameSetter[Dynamic[Filereadset], "Open"]
XlsWorksheed = "coefficient_FG_model";
XlsWorksheedThermo = "thermophysical data";
XlsWorksheedSpread = "parmeter_spread_model";
XlsWorksheedBond = "parameter_bond_strength_model";
XlsWorksheedFriction = "friction";
Panel[Column[{
InputField[Dynamic[XlsWorksheed]],
InputField[Dynamic[XlsWorksheedThermo]],
InputField[Dynamic[XlsWorksheedSpread]],
InputField[Dynamic[XlsWorksheedBond]],
InputField[Dynamic[XlsWorksheedFriction]]}]]
{Button["Import Koeffizienten:",
Importsheet = Import[Filereadset, {"Sheets", XlsWorksheed}];
ImportsheetThermo =
Import[Filereadset, {"Sheets", XlsWorksheedThermo}];
ImportsheetSpread =
Import[Filereadset, {"Sheets", XlsWorksheedSpread}];
Importsheetbond =
Import[Filereadset, {"Sheets", XlsWorksheedBond}];
Importsheetfriction =
Import[Filereadset, {"Sheets", XlsWorksheedFriction}]
], Dynamic[{XlsWorksheed, XlsWorksheedThermo, XlsWorksheedSpread,
XlsWorksheedBond, XlsWorksheedFriction}]
}
Importsheet // TableForm
Datenbank =
ToExpression[
Transpose[
Drop[Drop[Transpose[Importsheet[[5 ;; Length[Importsheet]]]],
1], -Length[
Transpose[Importsheet[[5 ;; Length[Importsheet]]]]] + 12]]];
ImportsheetThermo // TableForm
DatenbankTP =
ToExpression[
Transpose[
Drop[Drop[
Transpose[ImportsheetThermo[[5 ;; Length[ImportsheetThermo]]]],
1], -Length[
Transpose[
ImportsheetThermo[[5 ;; Length[ImportsheetThermo]]]]] + 14]]];
ThermoFkt[{a_, b_, c_, Mod_}] := Which[
Mod == 1, a + b (T + 273.15),
Mod == 2, a Exp[ b (T + 273.15)] + c,
Mod == 3, a + b (T + 273.15) + c (T + 273.15)^2
]
(*If[Mod<2,a +b T,a Exp[b T]+c,];*)
MaterialNummer := {2, 5};
(* {MatName,A,m1,m2,m3,m4,m5,m6,m7,m8,m9}=ToExpression[Datenbank[[\
MaterialNummer]]];*)
MaterialVek =
Table[Datenbank[[MaterialNummer[[zaehler]] ]], {zaehler, 1,
Length[MaterialNummer]}]
Style[TableForm[{{"Werkstoff", "A", "m1", "m2", "m3", "m4", "m5",
"m6", "m7", "m8", "m9"}, Transpose[MaterialVek]}],
FontFamily -> "Arial"]

初始条件：
EK = 0;
phivor = 0.01;
rA1 = 0.003;(* der Radius von Material A enthalten B an der \
z-Koordinate Z1 *)
rA0 = 0.0035;(* der Radius von Material A enthalten B an der \
z-Koordinate 0 *)
rB0 = 0.0031;(* der Radius von Material B für Fall a) an der \
z-Koordinate 0 *)
\[Mu]R = \[Mu]RAW =
0.1;(* der Reibwert zwischen Werkzeug und Material A *)
\[Mu]S = \
\[Mu]RAB = 0.4;(* der Reibwert zwischen Material A und B *)
\[Alpha] =
8 Degree;(* Ziehwinkel Ziehhol und äußerre Ziehwinkel von Material \
A *)
\[Alpha]B = 5 Degree;(* Ziehwinkel von Material B in Bereich 2 *)
vE = v1 =
0.01;(* Eintrittsgeschwindigkeit des Draht *)
cA = 835.02; \[Rho]A \
= 2764.35; \[Lambda]A = 216.95;(* spezifische Wärmekapazität, Dichte \
und Wärmeleitfähigkeit von Material A *)
TA0 = 293.;(* die Eintrittstemperatur von Material A in Bereich 1 *)(* \
Werkstoff A, Al99,5 *)
cB = 371.61; \[Rho]B = 9082.36; \[Lambda]B = \
423.14;(* spezifische Wärmekapazität, Dichte und Wärmeleitfähigkeit \
von Material B *)
TB0 = 293.;(* die Eintrittstemperatur von Material B in Bereich 1 *)(* \
Werkstoff B, SF-Cu *)
cW = 500.; \[Rho]W = 7860.; \[Lambda]W = 45.;(* \
spezifische Wärmekapazität, Dichte und Wärmeleitfähigkeit von \
Werkzeug *)
TW = 293.;(* die Eintrittstemperatur von Werkzeug in Bereich 1 *)(* \
Werkzeug,Stahl*)
TWA = 333.;(* die Temperatur des Werkzeugs *)
\
(*MaterialVek=Prepend[Transpose[MaterialVek],{{cA,cB},{\[Rho]A,\[Rho]\
B},{\[Lambda]A,\[Lambda]B}}]*)
RWA = Sqrt[cA \[Rho]A \[Lambda]A]/(
Sqrt[cA \[Rho]A \[Lambda]A] + Sqrt[cW \[Rho]W \[Lambda]W]);
RAB = Sqrt[cB \[Rho]B \[Lambda]B]/(
Sqrt[cB \[Rho]B \[Lambda]B] + Sqrt[cA \[Rho]A \[Lambda]A]);
hW = 0.02;(* Höhe der Aufheizzone im Werkzeug *)
\[Alpha]WA =
6000; (* Wärmeübergangskoeffizient zwischen Werkzeug und Material A \
*)
\[Sigma]fA[\[CurlyPhi]A_, \[CurlyPhi]pA_, TA_] :=
10^6 150.93 E^(-0.002647 TA) \[CurlyPhi]A^0.08938 E^(-0.0193/\
\[CurlyPhi]A) (1 + \[CurlyPhi]A)^(-0.00034 TA) E^(
0.03501 \[CurlyPhi]A) \[CurlyPhi]pA^0.0108;
\[Sigma]fB[\[CurlyPhi]B_, \[CurlyPhi]pB_, TB_] :=
10^6 425 E^(-0.000501 TB) \[CurlyPhi]B^0.24001 E^(-0.0000000394/\
\[CurlyPhi]B) \[CurlyPhi]pB^0.0701;
\[CurlyPhi]A00a =
N[EK Sqrt[3]/(
9 (rA0 + rB0)) (rA0 (2 Tan[\[Alpha]] + Tan[0]) +
rB0 (Tan[\[Alpha]] + 2 Tan[0]))] + phivor(* im Material A *);
\[CurlyPhi]B00b =
N[EK Sqrt[3]/(9 rB0) rB0 2 Tan[\[Alpha]B]] +
phivor(* im Material B *);
\[CurlyPhi]p1Aa[0] = (2 (rA0^2 - rB0^2)/(rA0^2 - rB0^2))/(
rA0^2 - rB0^2) vE (rA0 Tan[\[Alpha]] + rB0 Tan[0]);
\[CurlyPhi]p1Bb[0] = (rB0^2 vE)/rB0^3 Tan[\[Alpha]B];
\[Alpha]AB = (
2 \[Lambda]A \[Lambda]B)/(\[Lambda]B (rA0 -
rA1) + \[Lambda]A rB0);
MaterialVek =
Append[{cA, \[Rho]A, \[Lambda]A, cB, \[Rho]B, \[Lambda]B},
MaterialVek];

ZiehenBiplat[GeometrieVek_, ProcessVek_, MaterialVek_, NumVek_] :=
Block[{phiAvor, phiBvor, v1, vE, vz, TA0, TB0, alphaw, EK, rA, rA0,
rA1, Aa0, Ab0, Aa, Ab, zmax, \[Alpha], LGl, rB0, rB, NMat, phiA,
phiB, dotphivA, dotphivB, \[Micro]RAW, \[Micro]RAB, kF, cpA, rhoA,
lambdaA, cpB, rhoB, lambdaB,
KoeffizientenLoc, \[Lambda]AB, \[Alpha]AB, lsg, FktPlot,
FktPlotSigma, DGLVec,
zlI, \[CapitalDelta]q1EA, \[CapitalDelta]q1EB, dridz, driplus1dz,
q, qB, TA, TB, rA0loc, rA1loc, rBloc, vzAStart, vzA, v0A, v0B,
alphaloc},
{rA0, rA1, \[Alpha], LGl, rB0} = GeometrieVek;
{phiAvor, phiBvor, v1, TA0, TB0, alphaw, \[Micro]RAW, \[Micro]RAB} =
ProcessVek;
{cpA, rhoA, lambdaA, cpB, rhoB, lambdaB, KoeffizientenLoc} =
MaterialVek;
{EK} = NumVek;
(* Anzahl der Werkstoffe: hier Test mit 2 *)
NMat = Length[MaterialVek];
(* Ziehholgeoemtrie *)
zmax = Cot[\[Alpha]] (rA0 - rA1); (*
Ende des Ziehols und Übergang zur Glättlänge LGl *)
rA[z_, rA0loc_, rA1loc_, alphaloc_] :=
If[z < zmax, rA0loc - z Tan[alphaloc], rA1loc];
(* rA[z,rA0,rA1,\[Alpha]]:=rA0; *)
(* Berührungspunkt Material B mit Ziehhol *)
zlI = If[rB0 >= rA1,
z /. FindRoot[rA[z, rA0, rA1, \[Alpha]] - rB0 == 0, {z, 0}][[1]],
zmax + LGl];
(* Flächeninhalt Eintrittsquerschnitt *)
Aa0 = Pi (Power[rA0, 2] - Power[rB0, 2]);
Ab0 = Pi Power[rB0, 2];
(* momentaner Flächeninhalt *)
Aa[z_, rBloc_] :=
Pi (Power[rA[z, rA0, rA1, \[Alpha]], 2] - Power[rBloc, 2]);
Ab[z_, rBloc_] := Pi Power[rBloc, 2];
(* lokale Geschwindigkeit *)
(* Startlösung *)
vzAStart[z_, rB_] := v1 Aa[zmax, rB]/Aa[z, rB];
vzA[z_, rB_] := vzAStart[z, 0];
v0A := vzAStart[0, 0];
v0B := vzAStart[0, 0];
(* Umformgrad allgemein *)
phiA[Aa0_, Aa_, phiAvor_] := Log[Aa0/Aa[z, rBloc]] + phiAvor;
phiB[Ab0_, Ab_, phiBvor_] := Log[Ab0/Ab[z, rBloc]] + phiBvor;
(* Umformgeschwindigkeit allgmein *)
(* vz in [mm/s] ... \
Ziehgeschwindigkeit *)
dotphivA[ri_, riplus1_, dridz_, driplus1dz_, vz_] :=
2.0 Abs[vz (ri dridz + riplus1 driplus1dz) (rA0^2 - rB0^2)/
Power[(Power[ri, 2] - Power[riplus1, 2]), 2]];
dotphivB[ri_, riplus1_, dridz_, driplus1dz_, vz_] :=
Abs[vz (riplus1 driplus1dz) (rB0^2)/Power[(Power[riplus1, 2]), 2]];
(* Fließspannungsvektor FG-Ansatz*)
(* ProcessData: {phi,dotphi,T}*)
(* \
Koeffizienten: {MatName,A,m1,m2,m3,m4,m5,m6,m7,m8,m9} *)
(* Schicht: \
Nummer der Schicht *)
kF[ProcessData_, Koeffizienten_, Schicht_] :=
Koeffizienten[[Schicht,
2]] Exp[-Koeffizienten[[Schicht, 3]] ProcessData[[3]]] Power[
ProcessData[[1]] + 0.05, Koeffizienten[[Schicht, 4]]] Power[
ProcessData[[2]], Koeffizienten[[Schicht, 5]] ] Exp[
Koeffizienten[[Schicht, 6]]/(ProcessData[[1]] + 0.05)] Power[
1 + ProcessData[[1]] + 0.05,
Koeffizienten[[Schicht, 7]] ProcessData[[3]]] Power[
1 + ProcessData[[1]] + 0.05, Koeffizienten[[Schicht, 8]]] Exp[
Koeffizienten[[Schicht, 9]] (ProcessData[[1]] + 0.05)] Power[
ProcessData[[2]],
Koeffizienten[[Schicht, 10]] ProcessData[[3]]] Power[
ProcessData[[3]], Koeffizienten[[Schicht, 11]]];
\[Lambda]AB[z_] := ((rA[z, rA0, rA1, \[Alpha]] - rB[z])/2 \[Lambda]A +
rB[z]/2 \[Lambda]B)/(
rA[z, rA0,
rA1, \[Alpha]]/2);(* Wärmeleitfähigkeit mit z zwischen Material A \
und B *)
\[Alpha]AB = (
2 \[Lambda]A \[Lambda]B)/(\[Lambda]B (rA0 - rA1) + \[Lambda]A rB0);
(*Lösung im Bereich 1 für die Spannungs- und Temperaturgleichungen *)
If[kF[{phiAvor, 0.1, TA0 - 273.15}, KoeffizientenLoc, 1] <
kF[{phiBvor, 0.1, TB0 - 273.15}, KoeffizientenLoc, 2],
(* Schicht A weicher als Schicht B *)
rB[z_] := rB0;
dridz := Tan[\[Alpha]];
driplus1dz := 0;