adjustable number of differential equations

Hello,

I am trying to model non-radiative energy transfer between a single donor molecule towards a variable number of acceptor molecule that are on periphery of a circle. instead of having a fixed number of acceptor, I am trying to use the IntegrateODE command with a number that can be set by the user.

Function FRET_sphere(pw, tt, yw, dydt)
	Wave pw	// pw[0] = k1, pw[1] = k2, pw[2] = k3
	Variable tt	// time value at which to calculate derivatives
	Wave yw	// yw[0]-yw[3] containing concentrations of A,B,C,D
	Wave dydt	  // wave to receive dA/dt, dB/dt etc. (output)
	Variable /g na
	dydt[0] = -pw[0]*exp(-((tt)/0.03)^2)*yw[0]+pw[1]*yw[1] // pw[0]=donor excitation rate, pw[1]=1/(lifetime donor) 

	dydt[1] = pw[0]*exp(-((tt)/0.03)^2)*yw[0]-pw[1]*yw[1] - pw[na+2]*yw[1] //transfer toward na acceptors on a circle line with a total rate = pw[3]+...+pw[na+2]

variable i
for (i=0;i<(na);i=i+1)
	dydt[2+i] = pw[3+i]*yw[1]-pw[2]*yw[2+i] 
endfor

End

function EFRET_sphere(Rf,Racc,na)
variable Rf,Racc,na
variable i,j,r,aa,k
make /o/n=(na+4) pw
pw[0]=6500
pw[1]=0.25
pw[2]=0.33
make /o/d/n=20 FRET
SetScale/P x 1,1,"", FRET
FRET=0
Make/O/D/N=(200,na+2) ChemKin
SetScale/I x 0,50,"", ChemKin
Chemkin=0
ChemKin[0][0]= 1
for (k=1;k<2;k=k+1)
	for (j=0;j<(na);j=j+1)
		r=sqrt((Racc+k)^2+(Racc)^2-2*(Racc+k)*Racc*cos((pi/2)*(1-j/3)))	
		aa=aa+pw[1]*(Rf/r)^6
		pw[3+j]=pw[1]*(Rf/r)^6
	endfor
	pw[na+3]=aa
	IntegrateODE FRET_sphere, pw, ChemKin
	variable num=0
	wave out
	for (i=0;i<na;i=i+1)
		out[]=ChemKin[i+3][p]
		num=num+pw[2]*(sum(out))
	endfor
	out[]=ChemKin[1][p]
	variable den=num+pw[1]*sum(out)
	FRET[k]=num/den
endfor
end function

here is my code I suppose that IntegrateODE cannot handle a loop in the function defining the set of equations.

the second function EFRET_sphere is the part where I change the relative distance between the donor and the acceptors. 

Does anyone tried to use it in that way?

Best regards,

Pascal