gij={{1,0},{0,Sin[a]^2}}

gijk=ConstantArray[0,{2,2,2}];
Do[ Do[ Do[
gijk[[i,j,k]]=If[k==1,D[gij[[i,j]],a],D[gij[[i,j]],b]],{k,2}],{j,2}],{i,2}]

gIJ=Inverse[gij]
Glik=ConstantArray[0,{2,2,2}];
Do[ Do[ Do[
sum=Sum[gIJ[[l,j]] (gijk[[i,j,k]]+gijk[[j,k,i]]-gijk[[k,i,j]]),{j,2}];
Glik[[l,i,k]]=sum/2,
{k,2}],{l,2}],{i,2}]

Glik[[2,All,All]]                                                                                 
Out[44]= {{0, Cot[a]}, {Cot[a], 0}}

Glik[[1,All,All]]                                                                                 
Out[45]= {{0, 0}, {0, -(Cos[a] Sin[a])}}

{1/22}=-(Cos[a] Sin[a])
{2/12}={2/21}=Cot[a]

eqa=a''[t] -(Cos[a[t]] Sin[a[t]]) b'[t] b'[t]

eqb=b''[t] + 2 Cot[a[t]] b'[t] a'[t]

solution =NDSolve[{eqa==0,eqb==0,a[0]==Pi/2,b[0]==0,a'[0]==1,b'[0]==.5},{a,b},{t,0,5}]

ParametricPlot[Evaluate[{b[t],Pi/2-a[t]} /. solution],{t,0,5},PlotRange->{{0,2 Pi},{-Pi/2,Pi/2}}]

c3d[t_]=Evaluate[{Cos[b[t]] Sin[a[t]],Sin[b[t]] Sin[a[t]],Cos[a[t]]} /. solution]
ParametricPlot3D[c3d[t], {t, 0, 5}]

---
Glijk=ConstantArray[0,{2,2,2,2}];
Do[ Do[ Do[ Do [
Glijk[[l,i,j,k]]=Switch[k,1,D[Glik[[l,i,j]],a],2,D[Glik[[l,i,j]],b] ],{k,2}],{l,2}],{j,2}],{i,2}]

Rdabc=ConstantArray[0,{2,2,2,2}];

Do[ Do[ Do[ Do [
sum1=Sum[Glik[[j,a1,b1]] Glik[[d,c,j]]-Glik[[j,a1,c]] Glik[[d,b1,j]] ,{j,2}];
Rdabc[[d,a1,b1,c]]=FullSimplify[sum1+Glijk[[d,a1,b1,c]]-Glijk[[d,a1,c,b1]] ],{a1,2}],{b1,2}],{c,2}],{d,2}]

Rab=ConstantArray[0,{2,2}];
Do[ Do [
sum1=Sum[Rdabc[[j,a1,b1,j]],{j,2}];
Rab[[a1,b1]]=FullSimplify[sum1],{a1,2}],{b1,2}]

RAb=ConstantArray[0,{2,2}];
Do[ Do [
sum1=Sum[gIJ[[a1,j]] Rab[[j,b1]],{j,2}];
RAb[[a1,b1]]=FullSimplify[sum1],{a1,2}],{b1,2}]

R=Tr[RAb]

