ode 45 question, error message

I'm having trouble with my ode, it says

"Index exceeds the number of array elements (4)."

global p00
tic
a1 = 500;
a2 = 100;
a3 = 10^16.35
CA0 = 0.005; %FeS2
CB0 = 0.250;  %O2
CC0 = 0;     %Fe2
CD0 = 0.01;  %SO4
CE0 = 0.01; %Fe3
CAt = CA0 + CB0 + CC0 + CD0 + CE0;
CEt = CD0 + CC0 + CB0 + CE0
p00 = [a1;a2;a3];
y0(1) = CA0; % Concentration of initial A - [M]
y0(2) = CB0;
y0(3) = 0;
y0(4)= CD0;
yo(5) = CE0;

dt = 1;
t = (0:dt:60)';
tspan = [0 max(t)];
[T,y] = ode45('RK4odes',tspan,y0);
%
% Mass balance 
%
y(:,4) = CAt-y(:,1)-y(:,2)-y(:,3)-y(:,5)
y(:,5) = CEt-y(:,2)-y(:,3)-y(:,5)
%

toc


% seperate fxn
function RK4odes1= RK4odes(~,y)
global p00
a1 = p00(1);
a2 = p00(2);
a3 = p00(3);
RK4odes1(1,1) = -a1*(y(1)^0)*y(5)*(y(4)^-1)*10^-7 - a2*y(2)*(y(1)^0)*(y(4)^-1)*10^-7; % CA
RK4odes1(2,1) = -a2*(y(1)^0)*(y(3)^0)*(y(4)^-1)*10^-7 - a3*(y(3)^0)*y(2)*10^-7 % CB O
RK4odes1(3,1) = a2*(y(1)^0)*y(2)*(y(4)^-1)*10^-7 -a3*y(2)*(y(3)^0)*(y(4)^-1); % CC Fe2
RK4odes1(4,1) = a1*(y(1)^0)*y(5)*(y(4)^-1)*10^-7 +a2*y(2)*(y(1)^0)*(y(4)^-1)*10^-7;
RK4odes1(5,1) = -a1*(y(1)^0)*y(5)*(y(4)^-1)*10^-7 +a3*y(2)*(y(3)^0)*(y(4)^-1) 
end

Answer 1

[Star Strider]

This:

yo(5) = CE0;

should probably be:

y0(5) = CE0;

It likely does not help that ‘o’ and ‘0’ are next to each other on most keyboards!

[Estefania Garcia]

I didn't notice, thanks! I fixed it, but it still won't finish running

[Star Strider]

First, do not use global variables. Pass ‘p00’ as an additional argument:

function RK4odes1= RK4odes(t,y,p00)

Second, your system is ‘stiff’ with coefficients that vary significantly with respect to magnitude, so you must use a stiff solver:

[T,y] = ode15s(@(t,y)RK4odes(t,y,p00),tspan,y0);

Then, everything works!

On my Ryzen 9 3900 desktop:

Elapsed time is 0.050300 seconds.

.

[Estefania Garcia]

Thank you so much! This is an off topic question, but I'm new to coding with MATLAB and in general. How did you get better at it/what resources do you recommend?

[Star Strider]

My pleasure!

I have been involved in programming off and on since I took an undergraduate course in FORTRAN in 1968, and with MATLAB since 1993. Much of my current abilities I have developed while working with MATLAB Answers since 2012, with a few gaps to do other things. It is like everything else — I kept at it and learned from my own solutions and from reading others’ solutions. It just takes time, effort, and interest.

If my Answer helped you solve your problem, please Accept it!

.