%寫出一階常微分方程組
%ivpodefun.m
function dxdt =ivpodefun(t,x)
%初始化，22列1行矩陣，表示常微分方程組
dxdt = zeros(22,1);
%冷凝器微分方程
dxdt(1)=(17.4*1.763*x(2)/(1+0.763*x(2))-17.4*x(1))/100;
%精餾段微分方程組:從第1塊板到第9塊板
dxdt(2)=(17.4*1.763*x(3)/(1+0.763*x(3))-17.4*1.763*x(2)/(1+0.763*x(2))+13.05*(x(1)-x(2)))/25;
dxdt(3)=(17.4*1.763*x(4)/(1+0.763*x(4))-17.4*1.763*x(3)/(1+0.763*x(3))+13.05*(x(2)-x(3)))/25;
dxdt(4)=(17.4*1.763*x(5)/(1+0.763*x(5))-17.4*1.763*x(4)/(1+0.763*x(4))+13.05*(x(3)-x(4)))/25;
dxdt(5)=(17.4*1.763*x(6)/(1+0.763*x(6))-17.4*1.763*x(5)/(1+0.763*x(5))+13.05*(x(4)-x(5)))/25;
dxdt(6)=(17.4*1.763*x(7)/(1+0.763*x(7))-17.4*1.763*x(6)/(1+0.763*x(6))+13.05*(x(5)-x(6)))/25;
dxdt(7)=(17.4*1.763*x(8)/(1+0.763*x(8))-17.4*1.763*x(7)/(1+0.763*x(7))+13.05*(x(6)-x(7)))/25;
dxdt(8)=(17.4*1.763*x(9)/(1+0.763*x(9))-17.4*1.763*x(8)/(1+0.763*x(8))+13.05*(x(7)-x(8)))/25;
dxdt(9)=(17.4*1.763*x(10)/(1+0.763*x(10))-17.4*1.763*x(9)/(1+0.763*x(9))+13.05*(x(8)-x(9)))/25;
dxdt(10)=(17.4*1.763*x(11)/(1+0.763*x(11))-17.4*1.763*x(10)/(1+0.763*x(10))+13.05*(x(9)-x(10)))/25;
%進料板%精餾段微分方程組:從第10塊板
dxdt(11)=(17.4*(1.763*x(12)/(1+0.763*x(12))-1.763*x(11)/(1+0.763*x(11)))+13.05*x(11)-23.05*x(11)+5)/25;
%提餾段微分方程組:從第11塊板到第20塊板
dxdt(12)=(17.4*1.763*x(13)/(1+0.763*x(13))-17.4*1.763*x(12)/(1+0.763*x(12))+23.05*(x(11)-x(12)))/25;
dxdt(13)=(17.4*1.763*x(14)/(1+0.763*x(14))-17.4*1.763*x(13)/(1+0.763*x(13))+23.05*(x(12)-x(13)))/25;
dxdt(14)=(17.4*1.763*x(15)/(1+0.763*x(15))-17.4*1.763*x(14)/(1+0.763*x(14))+23.05*(x(13)-x(14)))/25;
dxdt(15)=(17.4*1.763*x(16)/(1+0.763*x(16))-17.4*1.763*x(15)/(1+0.763*x(15))+23.05*(x(14)-x(15)))/25;
dxdt(16)=(17.4*1.763*x(17)/(1+0.763*x(17))-17.4*1.763*x(16)/(1+0.763*x(16))+23.05*(x(15)-x(16)))/25;
dxdt(17)=(17.4*1.763*x(18)/(1+0.763*x(18))-17.4*1.763*x(17)/(1+0.763*x(17))+23.05*(x(16)-x(17)))/25;
dxdt(18)=(17.4*1.763*x(19)/(1+0.763*x(19))-17.4*1.763*x(18)/(1+0.763*x(18))+23.05*(x(17)-x(18)))/25;
dxdt(19)=(17.4*1.763*x(20)/(1+0.763*x(20))-17.4*1.763*x(19)/(1+0.763*x(19))+23.05*(x(18)-x(19)))/25;
dxdt(20)=(17.4*1.763*x(21)/(1+0.763*x(21))-17.4*1.763*x(20)/(1+0.763*x(20))+23.05*(x(19)-x(20)))/25;
dxdt(21)=(17.4*1.763*x(22)/(1+0.763*x(22))-17.4*1.763*x(21)/(1+0.763*x(21))+23.05*(x(20)-x(21)))/25;
%再沸器微分方程
dxdt(22)=(23.05*x(21)-17.4*1.763*x(22)/(1+0.763*x(22))-5.65*x(22))/200;