Main_AR_VandeVusse.m

% Construction of Attainable Region for van de Vusse Kinetics
% Author's Data: Housam BINOUS
% Department of Chemical Engineering
% National Institute of Applied Sciences and Technology
% Tunis, TUNISIA
% Email: [e-mail] 

global tau 

x0=[1 0];

[t,x]=ode45(@vandevusse,[0 10],x0);

Grad1=diff(x(:,2))./diff(x(:,1));

for i=1:1000
    r(i)=0.001*i;
    f(i)=interp1(x(1:220,1),Grad1(:,1),r(i))*(r(i)-1)-interp1(x(:,1),x(:,2),r(i));
end
f=f';
r=r';
z=interp1(f(100:980),r(100:980),-3e-7);

figure(1)

plot(x(:,1),x(:,2),'r')
axis([0 1 0 15e-5])
line([1 z],[0 interp1(x(:,1),x(:,2),z)])

j=1;
figure(2)
for i=1:10:10000
    tau=0.0001*i;
    options=optimset('display','off');
    X=fsolve(@cstr,[0.2 0.4],options);
    C(j,1)=X(1);
    C(j,2)=X(2);
    j=j+1;
end
j=1;
for i=1:100
    tau=i;
    options=optimset('display','off');
    X=fsolve(@cstr,[0.2 0.4],options);
    C2(j,1)=X(1);
    C2(j,2)=X(2);
    j=j+1;
end

Grad=diff(C(:,2))./diff(C(:,1));

for i=1:1000
    r(i)=0.001*i;
    f(i)=interp1(C(1:999,1),Grad(:,1),r(i))*(r(i)-1)-interp1(C(:,1),C(:,2),r(i));
end

z=interp1(f(100:980),r(100:980),-0.0004e-3);

plot(x(:,1),x(:,2),'r')
axis([0 1 0 15e-5])
hold on
plot(C(:,1),C(:,2),'g')
line([1 z],[0 interp1(C(:,1),C(:,2),z)])
plot(C2(:,1),C2(:,2),'g')

x0=[z interp1(C(:,1),C(:,2),z)];
[t,x]=ode45(@vandevusse,[0 10],x0);
plot(x(:,1),x(:,2),'m')