I have to fit to experimental data a model function given by a linear combination of functions y1(t) and y2(t)

        a1*y1(t) + a2*y2(t)

with a1 and a2 adjustable parameters.

The functions y1 and y2 are obtained by solving two coupled ODE's :

y1'(t)+y1(t)/t1-f(t)==0 with y1(0)==0

y2'(t)+y2(t)/t2-f(t)-a*y1(t)=0 with y2(0)==0

The function f(t) is known analytically.

Thus, there are 5 adjustable parameters a, a1, a2, t1, t2.

I found the following example in the Mathematica documentation, namely:

model[a_?NumberQ, b_?NumberQ, c_?NumberQ] :=  Module[{y, x},First[y /.NDSolve[{y''[x] + a y[x] == 0, y[0] == b, y'[0] == c}, 
 y, {x, 0, 10}]]]

nlm = NonlinearModelFit[data, model[a, b, c][x], {a, b, c}, x]

How can I generalize the given example to my case?