How to solve this ODE45 error?

clc
options=odeset('AbsTol', 1e-5, 'RelTol', 1e-8); 
sol2 = ode45(@sfun,[0, 100],[0.01 0.25  0 0 0 0 0 0 0 0 0 0 0]);
t=sol2.x;
p1=0.05*sol2.y(3,:)./sol2.y(1,:)
q1=0.3*sol2.y(5,:)./sol2.y(1,:)
s1=0.004*sol2.y(7,:)./sol2.y(1,:)
subplot(2,2,1), plot(sol2.x,p1,'b--', sol2.x,q1,'r--', sol2.x,s1,'k-.','linewidth',1)
xlabel('Time')
ylabel('B(n_a)_B(t,B)')
%axis([0 1.5 -4 0.5])
grid on
p2=0.05*sol2.y(4,:)./sol2.y(2,:)
q2=0.3*sol2.y(6,:)./sol2.y(2,:)
s2=0.004*sol2.y(8,:)./sol2.y(2,:)
subplot(2,2,2), plot(sol2.x,p2,'b--', sol2.x,q2,'r--', sol2.x,s2,'k-.','linewidth',1)
xlabel('Time')
ylabel('Bm_B(t,B)')
%axis([0 1.5 -0.05 0.01])
legend('B=\lambda','B=\beta','B=\eta')
grid on

%My function code is:

function v=sfun(t,z);


clc
lambda=0.03; beta=0.003; lambda0=0.009; mu=0.005;r=0.03; theta=0.4; omega=0.2; r0=0.01; p=0.3; eta=0.03; sigma=0.03;

v=[ sigma*(1-z(1))+ lambda*(z(2)/p+z(2))*(1-z(1))+ beta*(1-z(1))*z(1)-(lambda0+mu)*(1-z(1));
    eta*(1-z(1))+r*(1-theta*(z(1)/omega+z(1)))- r0*z(2);
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   -sigma*z(3)-(z(2)/p+z(2))*(1-z(1))-lambda*(z(2)/p+z(2))*z(3)+lambda*(1-z(1))*(p*z(4)/(p+z(2))^2)+ beta*(1-z(1))*z(3)- beta*z(1)*z(3)-(lambda0+mu)*z(3);
   -eta*z(3)+r*(1-theta(z(1)/omega+z(1)))*z(3)-r*z(1)*(theta*omega*z(3)/(omega+z(1))^2)- r0*z(4);
   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    -sigma*z(5)-lambda*(z(2)/p+z(2))*z(5)+lambda*(1-z(1))*(p*z(6)/(p+z(2))^2)+ beta*(1-z(1))*z(5)- (1-z(1))*z(1)-(lambda0+mu)*z(5);
    -eta*z(5)+r*(1-theta(z(1)/omega+z(1)))*z(5)-r*z(1)*(theta*omega*z(5)/(omega+z(1))^2)- r0*z(6);
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    -sigma*z(7)-lambda*(z(2)/p+z(2))*z(7)+lambda*(1-z(1))*(p*z(8)/(p+z(2))^2)+ beta*(1-z(1))*z(7)- beta*z(1)*z(7)-(lambda0+mu)*z(7);
    (1-z(1))- eta*z(7)+r*(1-theta(z(1)/omega+z(1)))*z(7)-r*z(1)*(theta*omega*z(7)/(omega+z(1))^2)- r0*z(8)];


%The program is giving following error.
Array indices must be positive integers or logical values.

Error in sfun (line 11)
   -eta*z(3)+r*(1-theta(z(1)/omega+z(1)))*z(3)-r*z(1)*(theta*omega*z(3)/(omega+z(1))^2)- r0*z(4);

Error in odearguments (line 92)
f0 = ode(t0,y0,args{:});   % ODE15I sets args{1} to yp0.

Error in ode45 (line 107)
  odearguments(odeIsFuncHandle,odeTreatAsMFile, solver_name, ode, tspan, y0, options, varargin);

Error in sansplot2 (line 3)
sol2 = ode45(@sfun,[0, 100],[0.049 0.052  0 0 0 0 0 0  0 0 0 0 0]);

Answer 1

[Torsten]

options=odeset('AbsTol', 1e-5, 'RelTol', 1e-8); 
sol2 = ode45(@sfun,[0, 100],[0.01 0.25  0 0 0 0 0 0]);
t=sol2.x;
p1=0.05*sol2.y(3,:)./sol2.y(1,:)
q1=0.3*sol2.y(5,:)./sol2.y(1,:)
s1=0.004*sol2.y(7,:)./sol2.y(1,:)
subplot(2,2,1), plot(sol2.x,p1,'b--', sol2.x,q1,'r--', sol2.x,s1,'k-.','linewidth',1)
xlabel('Time')
ylabel('B(n_a)_B(t,B)')
%axis([0 1.5 -4 0.5])
grid on
p2=0.05*sol2.y(4,:)./sol2.y(2,:)
q2=0.3*sol2.y(6,:)./sol2.y(2,:)
s2=0.004*sol2.y(8,:)./sol2.y(2,:)
subplot(2,2,2), plot(sol2.x,p2,'b--', sol2.x,q2,'r--', sol2.x,s2,'k-.','linewidth',1)
xlabel('Time')
ylabel('Bm_B(t,B)')
%axis([0 1.5 -0.05 0.01])
legend('B=\lambda','B=\beta','B=\eta')
grid on

%My function code is:
function v=sfun(t,z);
lambda=0.03; beta=0.003; lambda0=0.009; mu=0.005;r=0.03; theta=0.4; omega=0.2; r0=0.01; p=0.3; eta=0.03; sigma=0.03;
v=[ sigma*(1-z(1))+ lambda*(z(2)/p+z(2))*(1-z(1))+ beta*(1-z(1))*z(1)-(lambda0+mu)*(1-z(1));
    eta*(1-z(1))+r*(1-theta*(z(1)/omega+z(1)))- r0*z(2);
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    -sigma*z(3)-(z(2)/p+z(2))*(1-z(1))-lambda*(z(2)/p+z(2))*z(3)+lambda*(1-z(1))*(p*z(4)/(p+z(2))^2)+ beta*(1-z(1))*z(3)- beta*z(1)*z(3)-(lambda0+mu)*z(3);
    -eta*z(3)+r*(1-theta*(z(1)/omega+z(1)))*z(3)-r*z(1)*(theta*omega*z(3)/(omega+z(1))^2)- r0*z(4);
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    -sigma*z(5)-lambda*(z(2)/p+z(2))*z(5)+lambda*(1-z(1))*(p*z(6)/(p+z(2))^2)+ beta*(1-z(1))*z(5)- (1-z(1))*z(1)-(lambda0+mu)*z(5);
    -eta*z(5)+r*(1-theta*(z(1)/omega+z(1)))*z(5)-r*z(1)*(theta*omega*z(5)/(omega+z(1))^2)- r0*z(6);
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    -sigma*z(7)-lambda*(z(2)/p+z(2))*z(7)+lambda*(1-z(1))*(p*z(8)/(p+z(2))^2)+ beta*(1-z(1))*z(7)- beta*z(1)*z(7)-(lambda0+mu)*z(7);
    (1-z(1))- eta*z(7)+r*(1-theta*(z(1)/omega+z(1)))*z(7)-r*z(1)*(theta*omega*z(7)/(omega+z(1))^2)- r0*z(8)];
end    

[Akshita Bhardwaj]

Thank you for the reply. It really worked.

I did the same earlier too but I don't know why it did not work. Still thanks a lot :)