function dx_dt=column(t,x)

global alpha n_stage n_feed F_in z_in qf refluxi vapori ...
    md mb mt 


reflux=refluxi;
vapor=vapori;

lr=reflux;
ls=reflux+F_in*qf;

vs=vapor;
vr=vs+F_in*(1-qf);

dist=vr-reflux;
lbot=ls-vs;

for i=1:n_stage
    y(i)=(alpha*x(i))/(1.+(alpha-1.)*x(i));
end

dx_dt(1)=(1/md)*(vr*y(2)-(dist+reflux)*x(1));

for i=2:n_feed-1;
    dx_dt(i)=(1/mt)*(lr*x(i-1)+vr*y(i+1)-lr*x(i)-vr*y(i)); %cambio vs en y(i+1) por vr
end

dx_dt(n_feed)=(1/mt)*(lr*x(n_feed-1)+vs*y(n_feed+1)-ls*x(n_feed)-vr*y(n_feed)+F_in*z_in);

for i=n_feed+1:n_stage-1;
    dx_dt(i)=(1/mt)*(ls*x(i-1)+vs*y(i+1)-ls*x(i)-vs*y(i));
end

dx_dt(n_stage)=(1/mb)*(ls*x(n_stage-1)-lbot*x(n_stage)-vs*y(n_stage));

dx_dt=dx_dt';