// 3 large vats with volume V1, V2, and V2 are connected in a circle
// with vat 1 going to vat 2 going to vat 3 going to vat 1 at one volume
// unit per one time unit. Each vat is well stirred  
// Let x, y and z be the amount of salt in the water of vat 1, 2 and 3
//
// This can be modelled by the system
// dx/dt = -x/V1 + z/V3
// dy/dt = x/V1 - y/V2
// dz/dt = y/V2 - z/V3
// Initially x(0) = 900, y(0) = 10, z(0) = 90.
//
// Compare the z values for the first 20 time units
// V1 = 100, V2 = 75, V3 = 25
// V1 = 100, V2 = 50, V3 = 50
// V1 = 100, V2 = 25, V3 = 75
//
function wdot = f(t, w)
wdot = [ -w(1,1)/V1 + w(3,1)/V3; w(1,1)/V1 - w(2,1)/V2; w(2,1)/V2 - w(3,1)/V3];
endfunction
w0 = [900; 10; 90];
t0 = 0;
t = 0:1.0:300;
V1 = 100; V2 = 75; V3 = 25;
w25 = ode(w0,t0,t,f);
V1 = 100; V2 = 50; V3 = 50;
w50 = ode(w0,t0,t,f);
V1 = 100; V2 = 25; V3 = 75;
w75 = ode(w0,t0,t,f);
scf(1);
clf();
plot(t,w25(1,:),'ro-',t,w50(1,:),'bo-',t,w75(1,:),'mo-');
xtitle('red V3=25, blue V3=50, magenta V3=75: amount of salt in vat 1 vs time','time','salt');
scf(2);
clf();
plot(t,w25(2,:),'ro-',t,w50(2,:),'bo-',t,w75(2,:),'mo-');
xtitle('red V3=25, blue V3=50, magenta V3=75: amount of salt in vat 2 vs time','time','salt');
scf(3);
clf();
plot(t,w25(3,:),'ro-',t,w50(3,:),'bo-',t,w75(3,:),'mo-');
xtitle('red V3=25, blue V3=50, magenta V3=75: amount of salt in vat 3 vs time','time','salt');