adjustable number of differential equations

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