Quiet@Remove["Global`*"];
$HistoryLength=2;
SetDirectory[NotebookDirectory[]];

n=2*3*5;
d=1.5;
λ=1.7;
cp=840;
ρ=1600;
a=(λ/(ρ*cp));
Print["a = ",N[a,2],"  (m^2).(s^(-1))"];
Δx=(d/n);
Print["Δx = ",Δx,"  m"];
ClearAll[Δt];
Fo=((a*Δt)/((Δx)^2));
Δtmax=Δt/.Solve[Fo==0.5][[1]];
Print["Δtmax = ",N[Δtmax,2],"  s"];

To=20;
ρel[T_]:=1.5*(10^(-3))*(1+0.0005*(T-To));
R=0.1;
proud=120;
Skolma=Pi*(R^2);
obvod=2*Pi*R;
α=5.5;
ϵ=0.9;
σ=5.67*(10^(-8));
QV1[T_]:=((((proud/Skolma)))^2)*ρel[T];
QV2[T_]:=(obvod/Skolma)*(α*(T-To)+ϵ*σ*((((T+273.))^4)-(((To+273.))^4)));
qv[T_]:=((QV1[T]-QV2[T])/(ρ*cp));

Δt=200;
Tpoc=To;
Tpolepoc=Array[Tpoc&,n];
To1=To;
Ton=To;
α1=100;
αn=α1;
fT1[T2_]=T1/.Solve[((-λ)/(Δx))*(T2-T1)==α1*(To1-T1),T1][[1]];
fTn[Tnminus_]=Tn/.Solve[((-λ)/(Δx))*(Tn-Tnminus)==αn*(Tn-Ton),Tn][[1]];

ClearAll[krok];
krok[list_List]:=Module[{Tpom,odQV},
odQV=Δt*qv[#]&/@(list[[2;;-2]]);
Tpom=odQV+ListConvolve[{Fo,(1-2Fo),Fo},list];AppendTo[Tpom,fTn[Last[Tpom]]];PrependTo[Tpom,fT1[First[Tpom]]]];

hodina=3600;
tmax=1*hodina;
nsteps=Round[(tmax/(Δt))];
teppole=NestList[krok,Tpolepoc,nsteps];
maxT=Last[teppole]//Max;

plotik[list_List,ntime_Integer]:=Module[{nel=Length[list],listpom,i,tpopis},listpom=Table[{(i-1)*Δx,list[[i]]},{i,1,nel}];ListPlot[listpom,AxesLabel->{"x[m]","T[°C]"},PlotRange->{Automatic,{0,maxT+10}},tpopis=ntime*Δt;PlotLabel->"t="<>ToString[tpopis]<>" s",Joined->True,PlotStyle->{Red,Thick}]];

plots=Table[plotik[teppole[[i]],i-1],{i,1,Length[teppole]}];

(*Manipulate[plots[[i]],{i,1,Length[plots],1}]*)

rceUstStav={a*T''[x]+qv[T[x]]==0}//Simplify;
okr1={-λ*T'[0]==α1*(To-T[0])};
okr2={-λ*T'[d]==-α1*(To-T[d])};
(*resNDS=NDSolve[Union[rceUstStav,okr1,okr2],T,{x,0,d}]
Plot[T[x]/.resNDS,{x,0,d}]*)

(*DSolve[rceUstStav,T[x],x]*)

resStrelba[T0_?NumericQ]:=Module[{rce,res},
rce=Union[rceUstStav,{-λ*T'[0]==α1*(To-T0)},{T[0]==T0}];
res=NDSolve[rce,{T,T'},{x,0,d}][[1]];
{-λ*T'[d]+α1*(To-T[d]),T}/.res
]
resStrelba[20.1]//Timing

(*Plot[resStrelba[T0][[1]],{T0,20,20.1}]*)

resStrelba[19.8][[2]];
Plot[%[x],{x,0,d}]

der2[{Tvlevo_,Tzde_,Tvpravo_}]:=((Tvlevo+Tvpravo-2Tzde)/((Δx)^2));
rceUvnitr[trojice:{Tvlevo_,Tzde_,Tvpravo_}]:=a*der2[trojice]+qv[Tzde]==0;
teploty=(Subscript[T, #])&/@Range[n];
rcesIn=rceUvnitr/@Partition[teploty,3,1]//Simplify;
Length[rcesIn]+2==n

rceT1={((-λ)/(Δx))*((Subscript[T, 2])-(Subscript[T, 1]))-α1*(To1-(Subscript[T, 1]))==0};
rceTn={((-λ)/(Δx))*((Subscript[T, n])-(Subscript[T, n-1]))-α1*((Subscript[T, n])-Ton)==0};

rceSite=Union[rcesIn,rceT1,rceTn];

Solve[qv[Tustalene]==0]
Tust=Select[Tustalene/.%,Im[#]==0&]//Max

pocTeploty={#,Tust}&/@teploty;
FindRoot[rceSite,pocTeploty]