GNU Octave Script Help - ODE

I'm trying to solve the following ODE:

where R(T) is defined as:

This is my not so great attempt at using Octave:

1;

function xdot = f(t, T)
xdot = 987 * ( 0.0000696 * ( 1 + 0.0038 * ( T(t) - 25 ))) - ( 0.0168 * (T(t)-25 )) - (( 3.25 * 10 ^ (-13))) * ((T(t))^4 - (25^4));
endfunction

[x, istate, msg] = lsode( "f", 100, (t=linspace(0,3600,1000)'));

T_ref and T_infinity_sign are the same constant. Why isn't my code correct?

Solution

If you type

debug_on_error(1)

on your octave session, and then run your code, you will see that the "f" part is called as expected, but then it fails inside lsode with the following error:

error: T(100): out of bound 1 (dimensions are 1x1)
error: called from
    f at line 4 column 8

If you look at the documentation of lsode, it says it expects a function f whose first argument is a state vector x, and the second is a scalar, corresponding to time t at which that state vector occurs; f is expected to output the differential dx/dt at time t for state vector x.

From your code it seems that you are reversing the order of the arguments (and their meanings).

Therefore, when you passed T as a second argument to your function, lsode treats it like a scalar, so when you then try to evaluate T(t), it fails with an 'out_of_bounds_ error.

My advice would be, read the documentation of lsode and have a look at its examples carefully, and start playing with it to understand how it works.

_{PS. The debug_on_error directive launches the debugger if an error occurs during code execution. To exit the debugger type dbquit (or if you're using the GUI, click on the 'Quit Debugging Mode' button at the top right of the octave editor). If you don't know how to use the octave debugger, I recommend you spend some time to learn it, it is a very useful tool.}