% One dimensional diffusion
x=[];k=0;S=[];P=[];Sf=[];
m=5;%number of working cells
d=.3;v=.1;vf=.1;%process parameters

%The matrix
P=zeros(m+1,m+1);
for i=1:m;P(i,i+1)=d;P(i+1,i)=d+v;end;P(m,m+1)=0;P(m+1,m)=vf;
for i=1:m+1;P(i,i)=1-sum(P(:,i));end;%P(m+1,m)=1;P(m,m)=0;P(m-1,m)=0;
if min(min(P))<0;warning='Matrix is incorrect',break;end
%The initial state vector
S0=zeros(m+1,1);%S0(1)=1;
S(:,1)=S0;

N=299;
Sf=zeros(m+1,N);
for i=1:N
%The feed vector (select one, close by % others) 
%if i==1;Sf(1,i)=1;else Sf(1,i)=0;end;o=1;%instant injection
Sf(1,i)=1;o=3;%continuous constant feed
%Tf=60;Sf(1,i)=1+.2*sin(2*pi*i/Tf);o=5;%sinusoidal feed
%Evolution
S(:,i+1)=P*S(:,i)+Sf(:,i);
    end
Sd=diff(S(m+1,:));%Calculation of differential output - outflow
im=[1:N]*Sd';%mean residence time. Makes sense only for the feed as instant injection
figure(1);colormap(white);bar3(S(1:m,:));pbaspect([2 1 1]);axis([ .5 N+.5  .5 m+1+.5 0 1.1*max(max(S(1:m,:)))]);
figure(2);subplot(3,2,o);plot(Sf(1,:));hold off;subplot(3,2,o+1);plot(diff(S(m+1,:)));hold on

