There are sevral problems I see with your formulation. First, the
functions you defined must be of the following form:
f1[t_] := c Sin[t]
f2[t_] := d Cos[t]
f3[t_] := Sqrt[f1[t]^2 + f2[t]^2]
x[t_] := ArcCos[f1[t]/f2[t]]
y[t_] := g Sec[x[t]]
Note that when you specify f1,f2 as arguments to another function you
must use f1[t], f2[t] etc. Then when you evaluate y[t] you get
In[6]:=y[t]
Out[6]=(d*g*Cot[t])/c
Then I would craete an array that is a function of t
In[7]:=r[t_] = Array[a[t], {3, 3}]
Out[7]={{a[t][1, 1], a[t][1, 2], a[t][1, 3]},
{a[t][2, 1], a[t][2, 2], a[t][2, 3]},
{a[t][3, 1], a[t][3, 2], a[t][3, 3]}}
Next we define the recursive functions
In[11]:=s0 = {1, 0, 0};
sf[1] = r[1].s0
sf[n_] := r[n].sf[n - 1]
Out[12]={a[1][1, 1], a[1][2, 1], a[1][3, 1]}
Thus evaluating sf[1] we get
In[14]:=sf[1]
Out[14]={a[1][1, 1], a[1][2, 1], a[1][3, 1]}
and simlarly for any other time step. Finally we need to substitute in
for the values of a[t][i,j]. I suggest a set of replacement rules Thus
for your function a[t][1,1] we have
In[15]:=sf[1] /. a[t_][1, 1] -> Sin[b]^2*Cos[y[t]] + Cos[x[t]]^2
Out[15]=Cos[(d*g*Cot[1])/c]*Sin[b]^2 + (c^2*Tan[1]^2)/d^2,
a[1][2, 1], a[1][3, 1]}
If you create an array of replacement rules (call them myrules) for
the a[t][i,j] terms, then you can be applied the rules in the
definition of sf[n]. Finally to plot the values we can use ListPlot
ListPlot[Table[{i,N[sf[i]]}/.myrules,{i,1,10}]]
I hope I have understood your problem correctly, if not it was a happy
invention
Brian