global N1 N2 N3

N3=0;

N=fsolve('shooting',[0.001 0.001],optimset('TolFun',1e-13));

N1=N(1) 
N2=N(2)

[t,y]=ode45('MaxwellStefan',[0 1.0],[0.319 0.528]);

% plot the solution
figure(1);
hold on;
Axis([0 1.0 0 1]);
plot(t,y(:,1),'r')
plot(t,y(:,2),'b')
plot(t,1-y(:,1)-y(:,2),'g')

% write the data to ASCII file
% columns:
% 1: coordinate
% 2: species 1 mole fraction
% 3: species 2 mole fraction
% 4: species 3 mole fraction
data = [t y 1-y(:,1)-y(:,2)];
save('maxwell_stefan_reference.res', 'data', '-ASCII')