cp = 4186;
kS = 1000;
Δ1 = T3 - Tin;
Δ2 = T4 - T2;
Pz = 1000;
m1 = 1;
m0 = 0.1;
m2 = 3;
Tin = 4;
T4 = 30;
Tw = 50;
ϵ = 
        0.5*((Tw + 273)/(Tw - Tout));
rce1 = 
        m0*Tin + k*m2*T4 == 
          m0*T2 + k*m2*T3;
 (*dostala voda*) rce2 = 
        m0*cp*(T2 - Tin) == 
          kS*((Δ1 + Δ2)/2);
 (*miseni*) 
rce3 = 
        k*m2*T3 + (1 - k)*m2*T4 == 
          m2*T5;
  
rce4 = 
        Pel + m2*cp*(T5 - Tout) == 
          m1*cp*(Tw - T6);
rce5 = 
        m1*cp*(Tw - T6) == ϵ*
            Pel;
rce6 = 
        cp*(m1 - m0)*(Tw - T7) == 
          Pz;
rce7 = (m1 - m0)*T7 + m0*T2 == 
          m1*T6;
rcedohr[k_] = {rce1, rce2, rce3, rce4, rce5, rce6, 
          rce7};
ClearAll[res1];
res1[k_] := Solve[rcedohr[k]][[1]]

Plot[ϵ /. res1[k], {k, 0, 1}, 
        PlotStyle -> {Hue[0.63], Thickness[0.008]}, 
        AxesLabel -> {"k[1]", "ϵ[k]"}];

Plot[((Pz + cp*m0*(Tw - Tin))/Pel) /. res1[k], {k, 0, 1}, 
      PlotStyle -> {Hue[0.63], Thickness[0.008]}, 
      AxesLabel -> {"k[1]", "ϵTotal[k]"}]

T1 = 22;
T3 = -12;
c = 1010;
V = 7*12*3;
ρ = (341/(T1 + 273. ));
m1 = ((V*ρ)/(3*3600));
rce1 = 
        m1*(T1 - T2) == 
          m1*(T4 - T3);
rce2 = 
        c*m1*(T1 - T2) == 
          kS*((Δ1 + Δ2)/2);
Δ1 = T1 - T4;
Δ2 = T2 - T3;
rcedohr[kS_] = {rce1, rce2};
res[
          kS_] := Solve[
            rcedohr[kS]][[1]];
Pnutny = m1*c*(T1 - T3);
Prekuperace[kS_] := 
        m1*c*(T4 - T3) /. res[kS];
ucinnost[kS_] := 
        Prekuperace[kS]/Pnutny;
Plot[ucinnost[kS] /. res[kS], {kS, 0, 1000}, 
        PlotStyle -> {Hue[0.63], Thickness[0.008]}, 
        AxesLabel -> {"kS[W/K]", "η[1]"}];

c = 1010;
V = 7*12*3;
ρ = (341/(T1 + 273. ));
m = ((V*ρ)/(3*3600));
Pz = 100;
Tin = -12;
T2 = 22;

rce1 = m*c*(T1 - T2) == Pz;
rce2 = 
        k*m*(T2 - T6) == 
          m*(T5 - Tin);
rce3 = 
        m*c*(T5 - Tin) == 
          kS*(((Δ1 + Δ2)/2));
Δ1 = T6 - Tin;
Δ2 = T2 - T5;
rce4 = 
      k*T6 + (1 - k)*T2 == 
        T3;    
rce5 = 
      Pel + m*c*(T3 - T4) == 
        m*c*(T1 - T5);    
rce6 = 
      m*c*(T1 - T5) == ϵ*
          Pel;    
ϵ = 
      0.5*((T1 + 273)/(T1 - T4));

rcedohr[k_, kS_] = {rce1, rce2, rce3, rce4, rce5, 
          rce6};
ClearAll[res1];
res1[k_, 
          kS_] = Solve[
            rcedohr[k, kS], {T3, T6, T5, Pel, T4, 
              T1}][[1]];

Plot3D[Pel /. res1[k, kS], {k, 0.001, 1}, {kS, 0.001, 100}]