Index in position 2 exceeds array bounds (Index must not exceed 1) errors?

% Code to solve system of four simultaneous differential equation
clc
N0 = [9.66e17;0;0;0];   % Intital conditions of N1, N2, N3, N4 (9.66e17,0,0,0)
[T, N] = ode15s(@rate_eq,[0 126e-15], N0);  % Intital conditions of time from 0 to 126e-15
plot(T,N(:,4),T,N(:,3),T,N(:,2),T,N(:,1))   % Plotting of N1 N2 N3 N4 vs Time

function dN = rate_eq( t,N )
t = 126e-15;
W = [0.2716 0.4243 0.6110 0.8317 1.0862 1.3748 1.6972]*1e9;
PION = [0.7072 0.8874 1.0254 1.2375 1.4143 1.5911 1.7679]*1e11;
dN = zeros(4,7);    % Blank array to store N1 N2 N3 N4 from solution for each values of seven values of W and PION
Q = 10e10;      % Constant value
A = 10e8;       % Constant value
Ppd = 10e10;    % Constant value
% W and PION is not constant and its values are calculated in first part of
% code using FOR loop

for j = 1:length(W)
    dN(1,j) = (-W(j)*N(1,j)) + ((Q+A)*N(2,j)) + ((W(j)+Q+A)*N(3,j)) + 0;            % first equation 
    dN(2,j) = 0 + (-(Q+A)*N(2,j)) + ((Q+A)*N(3,j)) + 0;                             % second equation 
    dN(3,j) = (W(j)*N(1,j)) + 0 - (W(j)+(2*Q)+(2*A)+ Ppd + PION(j))*N(3,j) + 0;     % third equation 
    dN(4,j) = 0 + 0 + PION(j)*N(3,j) + 0;                                           % fourth equation 
end
end

%-----------end of code----------------

I want to store values of N1 N2 N3 N4 for each seven different values of W and PION.

I also want to plot N1 N2 N3 N4 vs time for each W and PION.

But I am getting the following errors:

Index in position 2 exceeds array bounds. Index must not exceed 1.

Error in test_soln_rate>rate_eq (line 24)

dN(1,j) = (-W(j)*N(1,j)) + ((Q+A)*N(2,j)) + ((W(j)+Q+A)*N(3,j)) + 0; % first equation

Error in odearguments (line 92)

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

Error in ode15s (line 148)

odearguments(odeIsFuncHandle, odeTreatAsMFile, solver_name, ode, tspan, y0, options, varargin);

Error in test_soln_rate (line 6)

[T, N] = ode15s(@rate_eq,[0 126e-15], N0); % Intital conditions of time from 0 to 126e-15

Answer 1

[Alan Stevens]

N0 is a 4x1 vector, but you call it with two indices in the for loop! Turn N(1,j) into N(1) and similarly for the other calls to N.

[JEESHAN ANSARI]

I have tried using N(1) and so on for other N values. But now I am getting the following errors:

Error using odearguments

RATE_EQ must return a column vector.

Error in ode15s (line 148)

odearguments(odeIsFuncHandle, odeTreatAsMFile, solver_name, ode, tspan, y0, options, varargin);

Error in test_soln_rate (line 6)

[T, N] = ode15s(@rate_eq,[0 126e-15], N0); % Intital conditions of time from 0 to 126e-15

[Alan Stevens]

Here's a quick and dirty way of solving the equations. I'll leave you to tidy it up!

% Code to solve system of four simultaneous differential equation
clc
N0 = [9.66e17;0;0;0];   % Intital conditions of N1, N2, N3, N4 (9.66e17,0,0,0)
dt = 126e-15/50;
tspan = 0:dt:126e-15;

for j = 1:7
    [T, N] = ode15s(@(t,N) rate_eq(t,N,j),tspan, N0);
    N1 = N(:,1); N2 = N(:,2); N3 = N(:,3); N4 = N(:,4);
    figure
    % Plotting of N1 N2 N3 N4 vs Time
    subplot(2,2,1)
    plot(T,N1),grid
    xlabel('T'),ylabel('N1') 
    subplot(2,2,2)
    plot(T,N2),grid
    xlabel('T'),ylabel('N2')
    subplot(2,2,3)
    plot(T,N3),grid
    xlabel('T'),ylabel('N3')
    subplot(2,2,4)
    plot(T,N4),grid
    xlabel('T'),ylabel('N4')
end

function dN = rate_eq(~,N,j )

    W = [0.2716 0.4243 0.6110 0.8317 1.0862 1.3748 1.6972]*1e9;
    PION = [0.7072 0.8874 1.0254 1.2375 1.4143 1.5911 1.7679]*1e11;
    w = W(j); p = PION(j);
    Q = 10e10;      % Constant value
    A = 10e8;       % Constant value
    Ppd = 10e10;    % Constant value
    N1 = N(1); N2 = N(2); N3 = N(3); N4 = N(4);
 
    dN1 = (-w*N1) + ((Q+A)*N2) + ((w+Q+A)*N3) + 0;            % first equation 
    dN2 = 0 + (-(Q+A)*N2) + ((Q+A)*N3) + 0;                   % second equation 
    dN3 = (w*N1) + 0 - (w+(2*Q)+(2*A)+ Ppd + p)*N3 + 0;    % third equation 
    dN4 = 0 + 0 + p*N3 + 0;                                % fourth equation 

    dN = [dN1; dN2; dN3; dN4];
end

[Alan Stevens]

@Torsten Thanks for the repair - I must avoid doing stuff when I'm tired!!

[JEESHAN ANSARI]

@Alan Stevens @Torsten

Thanks for helping.

[JEESHAN ANSARI]

I still have one problem remaining in the above code. I need seven columns of N1 N2 N3 N4 for each of seven values of W and PION.

But the code is giving only a single column of N1 N2 N3 N4 (I dont know for which value of W and PION).

[Alan Stevens]

Like this?

N0 = [9.66e17;0;0;0];   % Intital conditions of N1, N2, N3, N4 (9.66e17,0,0,0)
dt = 126e-15/50;
tspan = 0:dt:126e-15;
W = [0.2716 0.4243 0.6110 0.8317 1.0862 1.3748 1.6972]*1e9;
PION = [0.7072 0.8874 1.0254 1.2375 1.4143 1.5911 1.7679]*1e11;

for j = 1:7
    w = W(j); p = PION(j);
    [T, N] = ode15s(@(t,N) rate_eq(t,N,w,p),tspan, N0);
    N1(:,j) = N(:,1); N2(:,j) = N(:,2); N3(:,j) = N(:,3); N4(:,j) = N(:,4);
end
% Plotting of N1 N2 N3 N4 vs Time
figure
plot(T,N1),grid
xlabel('T'),ylabel('N1')
legend('j=1','j=2','j=3','j=4','j=5','j=6','j=7','Location','southwest')

figure
plot(T,N2),grid
xlabel('T'),ylabel('N2')
legend('j=1','j=2','j=3','j=4','j=5','j=6','j=7','Location','northwest')

figure
plot(T,N3),grid
xlabel('T'),ylabel('N3')
legend('j=1','j=2','j=3','j=4','j=5','j=6','j=7','Location','northwest')

figure
plot(T,N4),grid
xlabel('T'),ylabel('N4')
legend('j=1','j=2','j=3','j=4','j=5','j=6','j=7','Location','northwest')

function dN = rate_eq(~,N,w,p )

    Q = 10e10;      % Constant value
    A = 10e8;       % Constant value
    Ppd = 10e10;    % Constant value
    N1 = N(1); N2 = N(2); N3 = N(3); N4 = N(4);
 
    dN1 = (-w*N1) + ((Q+A)*N2) + ((w+Q+A)*N3) + 0;            % first equation 
    dN2 = 0 + (-(Q+A)*N2) + ((Q+A)*N3) + 0;                   % second equation 
    dN3 = (w*N1) + 0 - (w+(2*Q)+(2*A)+ Ppd + p)*N3 + 0;    % third equation 
    dN4 = 0 + 0 + p*N3 + 0;                                % fourth equation 

    dN = [dN1; dN2; dN3; dN4];
end

[JEESHAN ANSARI]

Thanks @Alan Stevens for updating the code for me again. Now I am getting N1 N2 N3 N4 arrays as desired.