indexing error in function code

Hi,

i am using the below codes but it showing the indexing error

function H = mass(t,y,P)
ph1_l1=P.ph1_l1;
ph1_l2=P.ph1_l2;
rm=P.rm;
rs=P.rs;
Cp=P.Cp;
ohm=P.ohm;
A1=P.A1;
p1=P.p1;
Ca=P.Ca;
rf=P.rf;
B=P.B;
Io=P.Io;
Ia=P.Ia;
Ib=P.Ib;
pc1=P.pc1;
Cs=P.Cs;
niu=P.niu;
fm=P.fm;
phi_lm=P.phi_lm;
I11=P.I11;
Zeeta=P.Zeeta;
C1=P.C1;
rM=P.rM;
N=P.N;
G1=P.G1;
Rl=P.Rl;
H = zeros(12,12);
H(1,1)=1;
H(2,2)=1;
H(2,8)=ph1_l1*rm;
H(2,10)=ph1_l2*rm;
H(3,3)=1;
H(4,4)=A*rm/rs;
H(4,8)=rm/rs;
H(5,5)=1;
H(6,6)=A*rm/rs;
H(6,10)=rm/rs;
H(7,7)=1;
H(8,8)=1;
H(9,9)=1;
H(10,10)=1;
H(11,11)=Cp*ohm;
H(12,12)=Cp*ohm;
end
function dydt = fu(t,y,P)
ph1_l1=P.ph1_l1;
ph1_l2=P.ph1_l2;
rm=P.rm;
rs=P.rs;
Cp=P.Cp;
ohm=P.ohm;
A1=P.A1;
p1=P.p1;
Ca=P.Ca;
rf=P.rf;
B=P.B;
Io=P.Io;
Ia=P.Ia;
Ib=P.Ib;
pc1=P.pc1;
Cs=P.Cs;
niu=P.niu;
fm=P.fm;
phi_lm=P.phi_lm;
I11=P.I11;
Zeeta=P.Zeeta;
C1=P.C1;
rM=P.rM;
N=P.N;
G1=P.G1;
Rl=P.Rl;
dydt = [y(2)
    -p1^4*y(1)-Ca*y(2)+ph1_l1*rm*rf^2*B*sin(rf*t)+ph1_l2*rm*rf.^2*B*sin(rf*t)+(rf.^2*B*sin(rf*t-rf*cos(rf*t))*Io)
    y(4)
    -A*y(3)*pc1^4-(Cs*pc1^4*A+Ca*A)*y(4)-Ca*y(8)*I11-y(11)*niu*(Ia-Ib)+fm*phi_lm+((rm/rs)*rf^2*B*sin(rf*t)-Ca*rf*Bcos(rf*t))*I11
    y(6)
    -A*y(5)*pc1^4-(Cs*pci^4*A+Ca*A)*y(6)-Ca*y(10)*I11-q(12)*niu*(Ia-Ib)+fm*phi_lm+((rm/rs)*rf^2*B*sin(rf*t)-Ca*rf*Bcos(rf*t))*I11
    q(8)
    f^2*ph1_l1*y(1)+Zeeta*ph1_l1*y(2)+2*C1*(rs/rM)*pc1^3*y(3)-f^2*y(7)-Zeeta*y(8)+rf^2*B*sin(rf*t)
    q(10)
    f^2*ph1_l2*y(1)+Zeeta*ph1_l2*y(2)+2*C1*(rs/rM)*pc1^3*y(5)-f^2*y(9)-Zeeta*y(10)+rf^2*B*sin(rf*t)
    -N*G1*y(4)-y(11)/Rl
    -N*G1*y(6)-y(12)/Rl];
end
P.ph1_l1=-0.012247303909990;
P.ph1_l2=0.635309628411612;
P.rm=1.509423275703383;
P.rs=1.960632785814804;
P.Cp=1.170700000000000e-07;
P.ohm=33.720923694042130;
P.A1=0.089447240667428;
P.p1=9.769092596233520;
P.Ca=0.001400000000000;
P.rf=0.5;
P.B=0.028248587570621;
P.Io=0.554857887911763;
P.Ia=-0.002864751824000;
P.Ib=0.952971893800000;
P.pc1=1.140600000000000;
P.Cs=0.010100000000000;
P.niu=9.520415341114659e-07;
P.fm=2;
P.phi_lm=0.503431090400000;
P.I11=0.2;
P.Zeeta=3.696152028525914e-04;
P.C1=0.839958085300000;
P.rM=4.241322704123775;
P.N=-0.006529965080390;
P.G1=1.363350182000000;
P.Rl=0.0002;
tspan=0:0.001:100;
y0=0;
opts = odeset('Mass',@(t,q) mass(t,q,P));
[t,q] = ode45(@(t,q) fu(t,q,P),tspan,y0,opts);

and the error that i am having

Index exceeds the number of array elements. Index must not exceed 1.

Error in fu (line 28)
dydt = [y(2)
    
Error in Msolve>@(t,q)fu(t,q,P) (line 30)
[t,q] = ode45(@(t,q) fu(t,q,P),tspan,y0,opts);

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 Msolve (line 30)
[t,q] = ode45(@(t,q) fu(t,q,P),tspan,y0,opts);

Answer 1

[Star Strider]

The initial subscript problem was the result of:

y0=0;

because there are 12 differential equations. There are several typographical errors that I did my best to correct, replacing ‘f’ with ‘fr’, ‘q’ with ‘y’, and ‘A’ with ‘A1’. Even with those changes and expanding ‘y0’, it will not run in the 55 seconds allowed here.

I encourage you to experiment with it —

P.ph1_l1=-0.012247303909990;
P.ph1_l2=0.635309628411612;
P.rm=1.509423275703383;
P.rs=1.960632785814804;
P.Cp=1.170700000000000e-07;
P.ohm=33.720923694042130;
P.A1=0.089447240667428;
P.p1=9.769092596233520;
P.Ca=0.001400000000000;
P.rf=0.5;
P.B=0.028248587570621;
P.Io=0.554857887911763;
P.Ia=-0.002864751824000;
P.Ib=0.952971893800000;
P.pc1=1.140600000000000;
P.Cs=0.010100000000000;
P.niu=9.520415341114659e-07;
P.fm=2;
P.phi_lm=0.503431090400000;
P.I11=0.2;
P.Zeeta=3.696152028525914e-04;
P.C1=0.839958085300000;
P.rM=4.241322704123775;
P.N=-0.006529965080390;
P.G1=1.363350182000000;
P.Rl=0.0002;
tspan=0:0.001:100;
% y0=0;
y0 = zeros(12,1);
opts = odeset('Mass',@(t,q) mass(t,q,P));
[t,q] = ode15s(@(t,q) fu(t,q,P),tspan,y0,opts);

figure
plot(t,q)
grid
legend(compose('y_%d', 1:12), 'Location','bestoutside')

function H = mass(t,y,P)
ph1_l1=P.ph1_l1;
ph1_l2=P.ph1_l2;
rm=P.rm;
rs=P.rs;
Cp=P.Cp;
ohm=P.ohm;
A1=P.A1;
p1=P.p1;
Ca=P.Ca;
rf=P.rf;
B=P.B;
Io=P.Io;
Ia=P.Ia;
Ib=P.Ib;
pc1=P.pc1;
Cs=P.Cs;
niu=P.niu;
fm=P.fm;
phi_lm=P.phi_lm;
I11=P.I11;
Zeeta=P.Zeeta;
C1=P.C1;
rM=P.rM;
N=P.N;
G1=P.G1;
Rl=P.Rl;
H = zeros(12,12);
H(1,1)=1;
H(2,2)=1;
H(2,8)=ph1_l1*rm;
H(2,10)=ph1_l2*rm;
H(3,3)=1;
H(4,4)=A1*rm/rs;
H(4,8)=rm/rs;
H(5,5)=1;
H(6,6)=A1*rm/rs;
H(6,10)=rm/rs;
H(7,7)=1;
H(8,8)=1;
H(9,9)=1;
H(10,10)=1;
H(11,11)=Cp*ohm;
H(12,12)=Cp*ohm;
end
function dydt = fu(t,y,P)
ph1_l1=P.ph1_l1;
ph1_l2=P.ph1_l2;
rm=P.rm;
rs=P.rs;
Cp=P.Cp;
ohm=P.ohm;
A=P.A1;
p1=P.p1;
Ca=P.Ca;
rf=P.rf;
B=P.B;
Io=P.Io;
Ia=P.Ia;
Ib=P.Ib;
pc1=P.pc1;
Cs=P.Cs;
niu=P.niu;
fm=P.fm;
phi_lm=P.phi_lm;
I11=P.I11;
Zeeta=P.Zeeta;
C1=P.C1;
rM=P.rM;
N=P.N;
G1=P.G1;
Rl=P.Rl;
dydt = [y(2)
    -p1^4*y(1)-Ca*y(2)+ph1_l1*rm*rf^2*B*sin(rf*t)+ph1_l2*rm*rf.^2*B*sin(rf*t)+(rf.^2*B*sin(rf*t-rf*cos(rf*t))*Io)
    y(4)
    -A*y(3)*pc1^4-(Cs*pc1^4*A+Ca*A)*y(4)-Ca*y(8)*I11-y(11)*niu*(Ia-Ib)+fm*phi_lm+((rm/rs)*rf^2*B*sin(rf*t)-Ca*rf*B*cos(rf*t))*I11
    y(6)
    -A*y(5)*pc1^4-(Cs*pc1^4*A+Ca*A)*y(6)-Ca*y(10)*I11-y(12)*niu*(Ia-Ib)+fm*phi_lm+((rm/rs)*rf^2*B*sin(rf*t)-Ca*rf*B*cos(rf*t))*I11
    y(8)
    rf^2*ph1_l1*y(1)+Zeeta*ph1_l1*y(2)+2*C1*(rs/rM)*pc1^3*y(3)-rf^2*y(7)-Zeeta*y(8)+rf^2*B*sin(rf*t)
    y(10)
    rf^2*ph1_l2*y(1)+Zeeta*ph1_l2*y(2)+2*C1*(rs/rM)*pc1^3*y(5)-rf^2*y(9)-Zeeta*y(10)+rf^2*B*sin(rf*t)
    -N*G1*y(4)-y(11)/Rl
    -N*G1*y(6)-y(12)/Rl];
end

EDIT — (20 Nov 2023 at 14:20)

Changed integrator from ode45 to ode15s.

.

[Imran]

Thank you sir, but as you mention here it's not working or not gives solution what can I do? Is ther any other method to solve multiple ODE

[Star Strider]

As always, my pleasure!

I just now noticed that because of the large magnitude range of the parameters, the system is ‘stiff’. I edited my answer to change the integrator to ode15s, then re-ran the corrected code, and it works!