Need help on using ODE45

I'm trying to solve the following ODE

dx(1)/dt = k1 (1 + k(2)x(3)) − x(1)x(2)^2 − k(3)x(1),

dx(2)/dt = x(1)x(2)^2 − x(2) + k(3)x(1),

dx(3)/dt = x(2) − k(4)x(3),

subject to the initial conditions x(0) = x0 and plot x1, x2 and x3 as a function of t using 10^4 points for 0 ≤ t ≤ T. It should take T and vectors x0 and k0 as inputs.

So far I have the following code:

function crosscatalator(x0,k,T)

sol = ode45(@(t,x) f(t,x,k),[0 T],x0);

t = linspace(0,T,10000); X = deval(sol,t); plot(t,X(1,:),t,X(2,:),t,X(3,:))

legend('x_1','x_2','x_3')

xlabel('t')

function f = f(~,x,k)

f = zeros(3,1);

f(1) = k(1)*(1+k(2)*x(3))-x(1)*x(2)^2-k(3)*x(1);

f(2) = x(1)*x(2)^2-x(2)+k(3)*x(1);

f(3) = x(2)-k(4)*x(3);

When I try with initial conditions x0=[1 1 1], k=[1 2 3 4], T=200. I get an error saying index is out of bounds. I've tried everything I can possibly do to fix this but I keep getting an error. What am I doing wrong?

[Michael Haderlein]

I get the plots, no error occurs. Do you maybe have a variable which is also called "crosscatalator"? If that's not the case, please tell us where exactly the error occurs.

[John Smith]

I don't have another variable called "crosscatalator". I have noted the error below. Thanks.

Answer 1

[Mischa Kim]

John, your code is good. Is the code above all contained within one file with name crosscatalator.m?

If so, try at the MATLAB command prompt

 clear; 
 x0 = [1 1 1]; k = [1 2 3 4]; T = 200;
 crosscatalator(x0,k,T)

and tell us if it is still not working.

[John Smith]

I'v tried it both all within one file and in separate files. They seem to give the same error.

This is what I get:

Attempted to access x(3); index out of bounds because numel(x)=2.

Error in ode45 (line 3)

f(1) = k(1)*(1+k(2)*x(3))-x(1)*x(2)^2-k(3)*x(1);

Error in crosscatalator (line 2)

sol = ode45(@(t,x) f(t,x,k),[0 T],x0);

I don't know why it seems to be working with both of you and not working for me.

[Mischa Kim]

Try the clear. I assume there is another x0 floating around somewhere that's a two-vector.

What do you get when you output x0 at the command window?

 >> x0

and within crosscatalator right after the function header (right before the ode45 call)? Use the debugger by setting a break point, see image below. Simply move your mouse over the dashed line to the right of the line number (3 here) and click. When you run the code execution will halt at the break point and you can hover the mouse over the variables for inspection.

[John Smith]

When I enter x0 at the command window I get "Undefined function or variable 'x0'.".

If I do it within crosscatalator before ODE45 call then it outputs the vector x0 = [1 1 1]. So it seems fine but I don't know what's wrong with it.

[John Smith]

Just noticed your edit. I tried that and it gives the x0 = [1 1 1], k=[1 2 3 4] and T=200. It's so strange why it's not working.

[Mischa Kim]

Can you attach the file(s) you are working with?

[John Smith]

This is it with all the code in the single file for convenience. I have tried the code both in the same file and separate files with no success. Thank you for helping.

[Mischa Kim]

One last (for now) change: re-name the ODE file. Of course, also in the function call, line 2.

 function f = myf(~,x,k)
 ...

[John Smith]

It's still the same with this change, doesn't work. I'm starting to think it's probably something wrong with matlab on the university computers...is that even possible?

I haven't tried testing the code using matlab on my laptop, so it might work as it worked for you, even so it seems strange that it's not working on university computers.