f1=inline('-3.9*m*t/(z^2)-1.2*m*a/(z^2)-14.5*m*b/(z^2)-1.34*m*s/(z^2)','m','t','a','b','s');
z=1+3.3a+13.6b+33.3f+215c
f2=inline('-3.9*m*t/(z^2)','m','t');
z=1+3.3a+13.6b+33.3f+215c
f3=inline('3.9*m*t/(z^2)-0.9*b/(z^2)-1.2*m*a/(z^2)','m','t','b','a')
z=1+3.3a+13.6b+33.3f+215c
f4=inline('2.3*b/(z^2)-1.34*m*s/(z^2)','b','m','s')
z=1+3.3a+13.6b+33.3f+215c
f5=inline('0.9*a/(z^2)-14.5*m*b/(z^2)-2.3*b/(z^2)','a','m','b')
z=1+3.3a+13.6b+33.3f+215c
f6=inline('1.2*a*m/(z^2)+14.5*m*b/(z^2)+1.34*m*f/(z^2)','a','m','b','s')
z=1+3.3a+13.6b+33.3f+215c
m(1)=3.8;
t(1)=6.2;
a(1)=0;
s(1)=0;
b(1)=0;
c(1)=0;
h=0.01;
y=0:h:20;
max=0;%max = Maximum for Selectivity of Para-xylene%
max2=0;%max2 = Maximum for Selectivity of Meta-xylene%
ObtimizedLength=1;
ObtimizedLength2=1;
for i=1:(length(y)-1)
k11 = hf1(m(i),t(i),a(i),b(i),s(i),c(i));
k12 = hf2(m(i),t(i),a(i),b(i),s(i),c(i));
k13 = hf3(m(i),t(i),a(i),b(i),s(i),c(i));
k14 = hf4(m(i),t(i),a(i),b(i),s(i),c(i));
k15 = hf5(m(i),t(i),a(i),b(i),s(i),c(i));
k16 = hf6(m(i),t(i),a(i),b(i),s(i),c(i));
k21 = h*f1(m(i)+k11/2,t(i)+k12/2,a(i)+k13/2,b(i)+k14/2,s(i)+k15/2,c(i)+k16/2);
k22 = h*f2(m(i)+k11/2,t(i)+k12/2,a(i)+k13/2,b(i)+k14/2,s(i)+k15/2,c(i)+k16/2);
k23 = h*f3(m(i)+k11/2,t(i)+k12/2,a(i)+k13/2,b(i)+k14/2,s(i)+k15/2,c(i)+k16/2);
k24 = h*f4(m(i)+k11/2,t(i)+k12/2,a(i)+k13/2,b(i)+k14/2,s(i)+k15/2,c(i)+k16/2);
k25 = h*f5(m(i)+k11/2,t(i)+k12/2,a(i)+k13/2,b(i)+k14/2,s(i)+k15/2,c(i)+k16/2);
k26 = h*f6(m(i)+k11/2,t(i)+k12/2,a(i)+k13/2,b(i)+k14/2,s(i)+k15/2,c(i)+k16/2);
k31 = h*f1(m(i)+k21/2,t(i)+k22/2,a(i)+k23/2,b(i)+k24/2,s(i)+k25/2,c(i)+k26/2);
k32 = h*f2(m(i)+k21/2,t(i)+k22/2,a(i)+k23/2,b(i)+k24/2,s(i)+k25/2,c(i)+k26/2);
k33 = h*f3(m(i)+k21/2,t(i)+k22/2,a(i)+k23/2,b(i)+k24/2,s(i)+k25/2,c(i)+k26/2);
k34 = h*f4(m(i)+k21/2,t(i)+k22/2,a(i)+k23/2,b(i)+k24/2,s(i)+k25/2,c(i)+k26/2);
k35 = h*f5(m(i)+k21/2,t(i)+k22/2,a(i)+k23/2,b(i)+k24/2,s(i)+k25/2,c(i)+k26/2);
k36 = h*f6(m(i)+k21/2,t(i)+k22/2,a(i)+k23/2,b(i)+k24/2,s(i)+k25/2,c(i)+k26/2);
k41 = h*f1(m(i)+k31,t(i)+k32,a(i)+k33,b(i)+k34,s(i)+k35,c(i)+k36);
k42 = h*f2(m(i)+k31,t(i)+k32,a(i)+k33,b(i)+k34,s(i)+k35,c(i)+k36);
k43 = h*f3(m(i)+k31,t(i)+k32,a(i)+k33,b(i)+k34,s(i)+k35,c(i)+k36);
k44 = h*f4(m(i)+k31,t(i)+k32,a(i)+k33,b(i)+k34,s(i)+k35,c(i)+k36);
k45 = h*f5(m(i)+k31,t(i)+k32,a(i)+k33,b(i)+k34,s(i)+k35,c(i)+k36);
k46 = h*f6(m(i)+k31,t(i)+k32,a(i)+k33,b(i)+k34,s(i)+k35,c(i)+k36);
m(i+1) = m(i) + (1/6)*(k11+2*k21+2*k31+k41);
t(i+1) = t(i) + (1/6)*(k12+2*k22+2*k32+k42);
a(i+1) = a(i) + (1/6)*(k13+2*k23+2*k33+k43);
b(i+1) = b(i) + (1/6)*(k14+2*k24+2*k34+k44);
s(i+1) = s(i) + (1/6)*(k15+2*k25+2*k35+k45);
c(i+1) = c(i) + (1/6)*(k16+2*k26+2*k36+k46);
R1(i+1) = f1(m(i+1),t(i+1),a(i+1),b(i+1),s(i+1),c(i+1));
R2(i+1) = f2(m(i+1),t(i+1),a(i+1),b(i+1),s(i+1),c(i+1));
R3(i+1) = f3(m(i+1),t(i+1),a(i+1),b(i+1),s(i+1),c(i+1));
R4(i+1) = f4(m(i+1),t(i+1),a(i+1),b(i+1),s(i+1),c(i+1));
R5(i+1) = f5(m(i+1),t(i+1),a(i+1),b(i+1),s(i+1),c(i+1));
R6(i+1) = f6(m(i+1),t(i+1),a(i+1),b(i+1),s(i+1),c(i+1));
spx(i)=a(i+1)/c(i+1);
if spx(i)>max
max=spx(i);
i=ObtimizedLength;
end
smx(i)=b(i+1)/c(i+1);
if smx(i)>max2
max2=smx(i);
i=ObtimizedLength2;
end
end
disp('max is = ');
disp(max);
disp('max2 is = ');
disp(max2);
disp('Obtimized Length is = ');
disp(ObtimizedLength);
disp('Obtimized Length2 is = ');
disp(ObtimizedLength2);
