%define a system
example_rxn = BioSystem();

%define rxn constants
example_rxn.AddConstant(Const('k1', 0.01));
example_rxn.AddConstant(Const('k2', 5));
example_rxn.AddConstant(Const('k3', 5));

%define the species in the rxn
A = Compositor('A', 100);
B = Compositor('B', 100);
C = Compositor('C', 0); %it's a good check to start
D = Compositor('D', 0); %products with zero
F = Compositor('F', 0); %as a sanity check

example_rxn.AddCompositor(A);
example_rxn.AddCompositor(B);
example_rxn.AddCompositor(C);
example_rxn.AddCompositor(D);
example_rxn.AddCompositor(F);

%create the parts
Part1 = Part('A+2B-->3C+D', [A B C D], ...      %elipsis lets you continue 
    [Rate('-const(k1)*val(A)*val(B)^2'), ...    %on the next line
    Rate('-2*const(k1)*val(A)*val(B)^2'), ...
    Rate('3*const(k1)*val(A)*val(B)^2'), ...
    Rate('const(k1)*val(A)*val(B)^2')]);
Part2 = Part('D+C-->F', [D C F], ...
    [Rate('-const(k2)*val(C)*val(D)'), ...
    Rate('-const(k2)*val(C)*val(D)'), ...
    Rate('const(k2)*val(C)*val(D)')]);
Part3 = Part('F-->D+C', [F D C], ...
    [Rate('const(k3)*val(F)'), ...
    Rate('const(k3)*val(F)'), ...
    Rate('-const(k3)*val(F)')]);

%add the parts to the system
example_rxn.AddPart(Part1);
example_rxn.AddPart(Part2);
example_rxn.AddPart(Part3);

%solve the differential equation
[T,Y] = example_rxn.run([ 0 0.25 ]); %determine run time interval

num_compositors = size(Y,2);
for i=1:num_compositors
	subplot(num_compositors, 1, i); %make a subplot for each compositor
	plot(T, Y(:,i));                %make the actual plot
	xlabel('time');ylabel(strcat('[',example_rxn.compositors(i).name,']'), 'Rotation', 0);
end