dest_cont.m

function varargout=dest_cont(varargin)
if nargin==0
Fplate=6;       %Plato de alimentacion
Np=20;          %Numero de platos incluido el hervidor y el condensador
F=1000;         %Alimentacion, kmol/h
q=1;            %Liquido saturado
xf=0.5;         %Composicion de la alimentacion 
R=2.25;         %Relacion de reflujo
alpha=2.46;     %Volatilidad relativa,  
MD=200;         %Masa en el condensador, kmol 
MR=400;         %Masa en el hervidor, kmol
M=50;           %Masa en los platos, kmol
V1=1575;        %Caudal vapor en la zona de agotamiento, kmol/h

else
    Fplate=varargin{1};
    Np=varargin{2};
    F=varargin{3};
    q=varargin{4};
    xf=varargin{5};
    R=varargin{6};
    alpha=2.46;
    MD=200;
    MR=400;
    M=50;
    V1=1575;
end

V=V1-(1-q)*F;  %Vapor en la zona de rectificacion, kmol/h
D=V/(R+1);     %Destilado, kmol/h
L=R*D;         %Liquido en la zona de rectificacion, kmol/h
L1=L+q*F;      %Liquido en la zona de agotamiento, kmol/h
W=L1-V1;       %Residuo, kmol/h

%Condicio0nes iniciales
X0=zeros(Np,1);


[t,x]=ode45(@ecuaciones,[0 10],X0,[],Fplate,Np,F,xf,alpha,MD,MR,M,V,L,L1,V1,W,D);

y=alpha*x./(1+(alpha-1)*x);

       
if nargin==0
subplot(2,2,1)
platos=1:Np;
mesh(platos,t,x)
xlabel('N Plato')
ylabel('Tiempo (h)')
zlabel('x')
title('Evolucin de la composicin en la fase lquida')
axis([1 Np 0 t(end) 0 1])
rotate3d on

subplot(2,2,2)
platos=1:Np;
mesh(platos,t,y)
xlabel('N Plato')
ylabel('Tiempo (h)')
zlabel('y')
title('Evolucin de la composicin en la fase vapor')
axis([1 Np 0 t(end) 0 1])
rotate3d on

subplot(2,2,3)
plot(t,x(1:end,1),'r')
xlabel('Tiempo (h)')
ylabel('xd')
title('Composicin en el destilado')
axis([0 t(end) 0 1])

subplot(2,2,4)
plot(t,x(1:end,Np),'r')
xlabel('Tiempo (h)')
ylabel('xw')
title('Composicin en cola')
axis([0 t(end) 0 max(x(1:end,Np)+2e-3)])
end

varargout{1}=t;
varargout{2}=x;
varargout{3}=y;
function dxdt=ecuaciones(t,x,Fplate,Np,F,xf,alpha,MD,MR,M,V,L,L1,V1,W,D)
   %Hervidor
   y(Np)=alpha*x(Np)/(1+(alpha-1)*x(Np));
   dxdt(Np,1)=(L1*x(Np-1)-V1*y(Np)-W*x(Np))/MR; 
 
   
for i=Np-1:-1:Fplate
    %Platos zona agotamiento
    y(i)=alpha*x(i)/(1+(alpha-1)*x(i));

    dxdt(i,1)=(L1*(x(i-1)-x(i))+V1*(y(i+1)-y(i)))/M;
   
end
    %Plato de alimentacion
    y(Fplate)=alpha*x(Fplate)/(1+(alpha-1)*x(Fplate));
    dxdt(Fplate)=(F*xf+L*x(Fplate-1)-L1*x(Fplate)+V*(y(Fplate+1)-y(Fplate)))/M;
   
for i=Fplate-1:-1:2
    y(i)=alpha*x(i)/(1+(alpha-1)*x(i));
    %Platos zona rectificacion
    dxdt(i,1)=(L*(x(i-1)-x(i))+V*(y(i+1)-y(i)))/M;
   
end
   i=1;

   %Composicion en el condensador
   y(i)=alpha*x(i)/(1+(alpha-1)*x(i));
   dxdt(1,1)=(V*y(i+1)-(L+D)*x(i))/MD;