solving odes

[Piotr]

Hello. I want to solve numerically the system of ODEs:

$y_0(t)=some\_given\_function,$

$y_n'(t)=y_n(t)a(t)+y_{n-1}(t)b(t)+y_n(t)x(t),  y_n(0)=y_n, n>0,$

where $a, b$ are given functions and $x$ is the solution of the 
auxiliary problem:

$x'(t)=f(t,x(t)), x(0)=x_0,$

and $n>0$ is some natural number, which can be sometimes huge.

I approximated solution of the auxiliary problem by lsode. Then I 
interpolated it by a piecewise linear function and solve the main 
equation for $n=1$ by lsode. For $n=2$ I tried to apply the same 
procedure with interpolation of $x$ and $y_1,$ but it failed. Could 
somebody suggest me a better (simpler) approach to this problem? I'm a 
newbie in Octave, so my idea is not sophisticated. I would be grateful 
for any help. Piotr.