File:Limit cycle of Brusselator.svg

clear;

//Parameters
a = 1.0;
b = 2.2;

// x=x(1), y=x(2)
function dx = BRU(t,x)
  dx(1) = a + x(1)^2 * x(2) - b * x(1) - x(1);
  dx(2) = b * x(1) - x(1)^2 * x(2); 
endfunction

// The time range and step for calculation
Tin = 0;
Tout = 40;
Tstep = 0.01;
t = Tin:Tstep:Tout;

// The initial values
x10 = 1;
y10 = 1.4;
X10 = [x10; y10];

X1 = ode(X10,0,t,BRU);

x20 = 1;
y20 = 2;
X20 = [x20; y20];

X2 = ode(X20,0,t,BRU);

x30 = 0.5;
y30 = 3;
X30 = [x30; y30];

X3 = ode(X30,0,t,BRU);

clf;

plot(X1(1,:), X1(2,:), 'r')
plot(X2(1,:), X2(2,:), 'b')
plot(X3(1,:), X3(2,:), 'g')

xlabel('x')
ylabel('y')