clear
echo on
%This script illustrates how a system of equations
%can be used to solve a chemical engineering problem.
%Consider a liquid-liquid extraction unit. In each stage
%of the unit the water and the solvent streams are
%assumed to be in equilibrium, such that y(i)=k*x(i)
%where k is the henry's law constant.  In addition,
%we have a mass balance on each stage which basically
%states that the amount of solute going into a stage
%equals to the amount going out - usually a good idea.
%
%We want to use this information to determine the
%concentration distribution in a four stage unit.
pause
%We let s be the molar flow rate of the solvent stream
%and we let w be the molar flow rate of the water
%stream.  k is the henry's law constant.  The
%solution depends solely on the ratio s*k/w.  We
%give this parameter the symbol chi.

chi = 1.5;
pause
%We set up the matrix:

n=4;
%This line of code creates the main diagonal of the
%matrix we wish to solve:
a=-(1+chi)*diag(ones(1,n),0)
pause

%This line creates the sub-diagonal of the matrix:
a=a+diag(ones(1,n-1),-1)
pause

%And this line creates the super diagonal:
a=a+chi*diag(ones(1,n-1),1)
pause

%We set up the right hand side:
b=zeros(n,1);
b(1)=-1
%where we have assumed that the inlet concentration
%of the solute in the water is unity, and that the
%inlet concentration of the solute in the solvent
%is zero.
pause

%We now solve for the concentration distribution
%using the \ command, and then we plot the result.
x=a\b;

figure(1)
plot(x,'o','MarkerSize',12)
set(gca,'FontSize',18)
axis([0.5 n+.5 0 1])
grid on
xlabel('stage number')
ylabel('concentration')
title('Extractor Concentration Distribution')
%Note that if only one vector is specified in the plot
%command, the vector (or matrix) is plotted against
%its index.  We also use the axis command to set the
%range of the axes in the plot to something which looks
%a bit better.

echo off