What is wrong with this code? Need help.

Main code:

clear all
close all
clc
M = [1 0 0
     0 1 0
     0 0 1];
x0=[-0.5; 0.2; 0.5];
options = odeset('Mass',M,'RelTol',1e-12,'AbsTol',[1e-14 1e-14 1e-14], 'Vectorized','on'); 
global t x y z dt alpha
dt=0.01;
for alpha=0.8800:0.0005:1.6000
alpha
clear n
clear m
[t,x]=ode15s(@equations,0:dt:500,x0);
n=length(x(:,1));
m=floor(n/2);
y=diff(x(m,n,1))/dt;   %% Error Message at this line: Index exceeds matrix dimensions
z=diff(y)/dt;
k=1;
clear aa
for i=m:n
t0=t(i); %Comp. local max. pts for m<t<n
option=optimset('display','off');
zer=fsolve(@differ,t0,option);
if interp1(t(m:n-2),z,zer)>0
aa(k)=interp1(t(m:n),x(m:n,1),zer);
k=k+1;
end
end
kmax=k-1;
h=plot(alpha.*ones(1,kmax),aa,'r.');
hold on
set(h,'MarkerSize',0.1);
end

Function equations.m:

function xdot=equations(t,x)
global alpha
a=0.0005; b=0.01; eps=0.01; beta=-1; R=0.3;
xdot(1) = (-x(2) + alpha*x(1)^2 + beta*x(1)^3)/eps;
xdot(2) = x(1) - x(3) - R*x(2);
xdot(3) = a - b*x(2);
xdot = xdot';
end

Function differ.m:

function f=differ(a)
global t x dt
s=diff(x(m:n,1))/dt;
f=interp1(t(m:n-1),s,aa);
end

The author of this program plot the attached picture. But when I tried cannot get that plot. Please help!!!!

Answer 1

[KALYAN ACHARJYA]

y=diff(x(m,n,1))/dt;   %% Error Message at this line: Index exceeds matrix dimensions

Yes, see the reason:

>> m
m =
       25000
>> n
n =
       50001
>> whos x
  Name          Size              Bytes  Class     Attributes
  x         50001x3               1200024  double    global    

The size of x is 5000x1 and you are trying to acess of x as follows (In first iterations)

x(25000,50001,1) 
%..^.........^ no x values 

Thats why it shows the error index exceeds matrix dimensions. That means you trying to acess the indices more than it have.

[ZAINULABADIN ZAFAR]

Respected Sir,

Can you fix it and plot the graph as attached in the pdf.

Thanks

[Walter Roberson]

Considering your function differ() uses

s=diff(x(m:n,1))/dt;

then I speculate that your earlier line

y=diff(x(m,n,1))/dt;

should be

y=diff(x(m:n,1))/dt;

[ZAINULABADIN ZAFAR]

Respected Sir,

I did tried for that also. It gives error. Can you help and plot. Plot is given in pdf file.

Thanks