syms Fj Vj V Area h I R Cpj Cp rhoj rho dH Tj Tji T Ck
f1 = Fj/Vj * (Tji - Tj) + (h*Area*(T-Tj))/(Vj*rhoj*Cpj);
f2 = (I * dH * exp(-Ck/(273 + T)) + I^2 * R - h * Area * (T-Tj))/(rho * Cp * V);
F = [f1;f2];
X = [T;Tj];
U = [Tji;Fj;I];
G = X;
A = sym(zeros(2,2));
B = sym(zeros(2,3));
C = sym(zeros(2,2));
D = sym(zeros(2,3));
for i = 1:2
    for j = 1:2
        A(i,j) = diff(F(i),X(j));
    end
end
for i = 1:2
    for j = 1:3
        B(i,j) = diff(F(i),U(j));
    end
end
for i = 1:2
    for j = 1:2
        C(i,j) = diff(G(i),X(j));
    end
end
for i = 1:2
    for j = 1:3
        D(i,j) = diff(G(i),U(j));
    end
end
A_subs = subs(A,[Fj Vj V Area h I R Cpj Cp rhoj rho dH Tj Tji T Ck],[108 30 50 30 0.00125 3 0.05 0.00097694 0.00027778 1.07 1.8 15000 30.47 25 46.93 4000]);
B_subs = subs(B,[Fj Vj V Area h I R Cpj Cp rhoj rho dH Tj Tji T Ck],[108 30 50 30 0.00125 3 0.05 0.00097694 0.00027778 1.07 1.8 15000 30.47 25 46.93 4000]);
C_subs = subs(C);
D_subs = subs(D);
A_dec = eval(A_subs);
B_dec = eval(B_subs);
C_dec = eval(C_subs);
D_dec = eval(D_subs);
ssmodel = ss(A_dec,B_dec,C_dec,D_dec);
tfmodel = tf(ssmodel);
gp23 = tfmodel(1,3);
gp23num = gp23.numerator{1};
gp23denom = gp23.denominator{1};