Autonomous linear differential equation with varying coefficients

Hi everybody, I need your kind suggestion about the following issue. I'm concerned in solving a differential equation in the function y(t). With y' and y'' I refer to the first and second time-derivative respectively. The equation is: y'' + a(y)y' + b(y)y = 0 y(0) = 2.1; y'(0) = 0; a, b are known per each value of the vector ips = [1.5,2], i.e. the configuration space which y belongs to. I tried the following strategy: I defined the vector tspan and called the ode45 solver;

tspan = 0:0.001:5;

ic = [2.1000 0];

opt=odeset('MaxStep',0.001);

[t,x] = ode45(@(t,x) myfun(t,x,ips,a,b), tspan, ic, opt);

the function myfun is:

function dxdt = myfun(t,x,ips,a,b)

dxdt = zeros(2,1) ;

a = interp1( ips, a, x(1),'linear' );

b = interp1( ips, b, x(1),'linear' );

dxdt(1,1) = x(2);

dxdt(2,1) = - b.*x(1) - a.*x(2) ;

This scheme does not work, since from some t on, solution is NaN.

Can anybody help me?

Thanks in advance

[Torsten]

And you are sure that interpolation makes sense, i.e. min(ips)<=x(1)<=max(ips) for all times t ?

Because x(1) at time t=0 is > max(ips) ...

Best wishes

Torsten.

[Giancarlo Santamato]

Hi Torsten, Thanks for answering. Actually I was wrong when I wrote ips, because ips = [1.5,2.195]. Your observation is correct: I'm not sure the interpolation makes sense; I followed the physical interpretation: ips is in fact the stroke of a slider, whose coordinate is y(t), while a and b are the y-dependent damping coefficient and square of the natural pulsation of a single degree of freedom kinematism. The equilibrium position I'm interested in is at 2.166, so if everything is correctly written, y(t) cannot overcome ips boundaries. I want to simulate the free response to a perturbation so as to test the possibility to model such a not-constant-parameters mechanism. Have you got any idea? Thank you so much