<< Calculus`VectorAnalysis`
SetCoordinates[Cartesian[x, y, z]];
rce = Div[Grad[T[x]]] == 0

resT = DSolve[{rce, T[0] == T1, T[d] == T2}, T[x], 
            x][[1]];
teplota = T[x] /. resT
qvect = -λ*Grad[teplota]
qx = qvect[[1]]
qx /. {d ->  .60, λ -> 0.8, T1 -> 15, 
        T2 -> 0}

dosabPha = Solve[{a*10 + b == (1/35), 
            a*50 + b == (1. /10)}][[1]]

dosabMB = Solve[{a*20 + b == (1/72), 
            a*75 + b == (1. /30)}][[1]]

ρ = 1000;
cp = 4186;
Ecko = (230/0.1);
rce = ρ*cp*T'[t] == (a*T[t] + b)*(Ecko^2) /. 
          dosabPha;
poc = T[0] == 10;
resT = DSolve[{rce, poc}, T[t], 
            t][[1]];
teplota = T[t] /. resT;
tmax = 
        t /. FindRoot[
            teplota == 100, {t, 100}];
Plot[teplota, {t, 0, tmax}, 
        AxesLabel -> {"t[s]", "T[°C]"}];
Print["tmax= ", (tmax/60), "  min"]

msteckou = 0.5;
cp = 4186;
T1 = 4;
T2 = 44;
ρ = 1000;
Print["P= ", ((msteckou*cp*(T2 - T1))/1000.) , " kW"]

Vsteckou = (msteckou/ρ);
y = 0.15;
d = 0.003;
Ecko = (230/d);
vx = (Vsteckou/(y*d));
rce = ρ*cp*vx*T'[x] == (a*T[x] + b)*(Ecko^2) /. 
          dosabPha;
teplota = 
        T[x] /. DSolve[{rce, T[0] == T1}, T[x], 
              x][[1]];
delka = 
        x /. FindRoot[
            teplota == T2, {x, 0.5}];
Plot[teplota, {x, 0, delka}, 
        AxesLabel -> {"x[m]", "T[°C]"}, 
        PlotStyle -> {Hue[0.63], Thickness[0.008]}];

Vsteckou = (msteckou/ρ);
y = 0.15;
d = 0.003;
Ecko = (230/d);
vx = (Vsteckou/(y*d));
rce = ρ*cp*vx*T'[x] == (a*T[x] + b)*(Ecko^2) /. 
          dosabPha;
teplota = 
        T[x] /. DSolve[{rce, T[0] == T1}, T[x], 
              x][[1]];
delka = 
        x /. FindRoot[
            teplota == T2, {x, 0.5}];
Plot[teplota, {x, 0, delka}, 
        AxesLabel -> {"x[m]", "T[°C]"}, 
        PlotStyle -> {Hue[0.63], Thickness[0.008]}];

a = 0.5;
b = 0.5;
c = 0.5;
cm = 0.01;
d1 = 3  cm;
Tpec = 1200;
S1 = 
        2*((a + d1)*(b + d1) + (c + d1)*(b + d1) + (a + 
                    d1)*(c + d1));
R1 = (d1/(λ1*S1));
datλ1 = {{1100, 0.4}, {1000, 0.36}} // 
          Union;
  
λ1 = 
        Fit[datλ1, {1, T}, T] /. 
          T -> 1200;
\(Print["λ1= ", λ1, \*"" W.(m^(-1)).(K^(-1))""]\)

S2 = 
        2*((a + 2  d1 + d2)*(b + 2  d1 + d2) + (c + 2  d1 + 
                    d2)*(b + 2  d1 + d2) + (a + 2  d1 + d2)*(c + 
                    2  d1 + d2));
\(\(Print["S2= ", 
        S2, \*"" \!\(m\^2\)""];\)\)
λ2 = 0.162;
R2 = (d2/(λ2*S2));

S3 = 
        2*((a + 2  d1 + 2  d2)*(b + 2  d1 + 2  d2) + (c + 2  d1 + 
                    2  d2)*(b + 2  d1 + 2  d2) + (a + 2  d1 + 
                    2  d2)*(c + 2  d1 + 
                    2  d2));
\(\(Print["S3= ", S3, \*"" \!\(m\^2\)""];\)\)

α = 5;
ϵ = 0.9;
σ = 5.67*(10^(-8));
Tor = 15;
Tok = 25;
Qk = ((Tpl - Tok)/((1/(S3*α))));
Qr = ϵ*σ*
          S3*((((Tpl + 273))^4) - (((Tor + 273))^4));

Tpl = 50;
rce1 = Q == Qk + Qr;
rce2 = Q == ((Tpec - Tpl)/(R1 + R2));
Solve[{rce1, rce2}][[3]]