MATLAB: Bifurcation diagram for a non autonomous system.

graphodeode45

I am trying to draw a bifurcation diagram for a non-autonomous system. The code:
function LKEbif
global b P theta c a e q d alpha1 Ke omega b=1.2;P=0.025;theta= 0.03; c= 0.81;a= 0.25; e=0.8; q=0.0038; d=0.25; alpha1=0.5;omega= 2*pi/6;
option=odeset('AbsTol',1e-10,'RelTol',1e-10);
limit=[1400:1:1500];
ymin=[];
ymax=[];
for Ke=0:0.01:3;
[~,U]=ode45(@Onetwo,[0:1:1500],[0.5;0.25;1.5]);
u2=U(limit,2);
ymin=[ymin,min(u2)];
ymax=[ymax,max(u2)];
end
X=0:0.01:3;
figure
plot(X, ymin,'g.',X, ymax,'k.');
legend('grazer min','grazer max'),
hold on
xlim([0 3]);
title('Bifurcation diagram for LKE model','FontSize',12)
xlabel('P, the total phosphorous');
ylabel('producer, grazer','Rotation',90);
function [Ydot] = Onetwo(t,y)
global b P theta c a e q d alpha1 Ke omega
Ydot(1) = b*y(1)*(1-y(1)/(min(y(3),(P-theta*y(2))/q)))-c*y(1)*y(2)/(a+y(1));
Ydot(2) =e*min(1,((P-theta*y(2))/y(1))/theta)*c*y(1)*y(2)/(a+y(1))-d*y(2);
Ydot(3)=(alpha1*(Ke-y(3)))+(Ke*omega*cos(omega*t)/2);
Ydot=Ydot';
The system gives me the following error:
Index exceeds matrix dimensions.
Error in LKEbif (line 24)
u2=U(limit,2);

Best Answer

Hi!
The reason for the error "Index exceeds matrix dimensions." in your piece of code is because U is an array of size 51 x 3 but you are trying to use 'limit' which range from 1400 to 1500 to index into U.
U(a,b) % will be valid only if a<52 and b<4 where a and b are positive integers OR Logicals.
In cases of trivial issues like this, I would suggest you to try debugging using breakpoints in MATLAB.
Hope this helps.
Regards,
Jyotish