global r s x1 I 

r = 0.001; 
s = 4; 
x1 = -(1+sqrt(5))/2; 
I = 2; 

[T,Y] = ode23('func_bursts',[0,1000],[-1.4485,-9.2475,2.1105]); 

figure(1) 
subplot(3,1,1) 
plot(T,Y(:,1)); 
xlabel('t'); 
ylabel('x'); 

subplot(3,1,2); 
plot(T,Y(:,2)); 
xlabel('t'); 
ylabel('y'); 

subplot(3,1,3); 
plot(T,Y(:,3)); 
xlabel('t'); 
ylabel('z'); 
axis([0 1000 1.5 2.5]) 
