Matlab ODE45: solving with coefficients that are functions of time

[hasidim]

Hello all,

I am new to the ODE solvers in Matlab and am trying to learn:

First, I am solving a 2nd order ODE to determine x(t), x'(t), and x''(t). No problem. Then, I am using these solutions to calculate two coefficients (that are functions of time) that are used in a second, second order ODE. The output solution is a vector of NaN's.

To explain... the first ODE script:

% Initial Conditions:
x10 = x0;  %x0 is some initial value
x20 = 0;  
tspan = linspace(0,tmax+dt, tmax/dt+1);

% Solve ODE
[t,x] = ode23('ode', tspan, [x10,x20]);

and first ODE function (where c1 and c2 are arbitrary, constant coefficients):

function xp = ode(t,x)
xp = zeros(2,1);
xp(1) = x(2); 
xp(2) = c2*x(1)-c1*x(2)

So now I have solutions for x(t), x'(t), x"(t) (that I name x, xp, xp2 respectively).

Ok, so no problem until I get to the second ODE. Following the same approach:

Script:

% Initial Conditions:
a10 = a0;  %a0 is some initial value
a20 = 0;  
tspan = linspace(0,tmax+dt, tmax/dt+1);

% Solve ODE
[t,a] = ode23('ode2', tspan, [a10,a20]);

ODE function:

function ap = ode2(t,a)

K1 = xp2./x;
K2 = xp1./x; % the length of K1 and K2 is the length of x from solution above

ap = zeros(2,1);
ap(1) = a(2); 
ap(2) = K2*a(1)+ K1*a(2)

Using the above ODE function, I will get an indicie mismatch error. If I try to apply indices to K1 and K2 [(floor(t/dt)+1)], the output will all be NaN's.

It seems like there should be a simple answer as to how to solve a second order ODE that has coefficients (that are functions of time) that I have already solved for. Thanks in advance. (Hopefully my explanation is clear).

 

[hasidim]

Hi All,

I found my issue. It was with the coefficient calculation, not with the ODE solver.

For what its worth, the coefficients in the second ODE (the coefficients that are functions of time) do need to be indexed as a function of t. So, something like this:

function ap = ode2(t,a)

K1 = xp2./x;
K2 = xp1./x; % the length of K1 and K2 is the length of x from solution above

ap = zeros(2,1);
ap(1) = a(2); 
ap(2) = K2(t/dt+1)*a(1)+ K1(t2/dt+1)*a(2)