Michael,

ParametricNDSolve already does much of what you want. You first need to clear out values you previously set

Clear[a, b, i, u, o]

Solve your problem:

sol = ParametricNDSolve[{x'[t] == (Sin[y[t]] - b*x[t])*
     Cos[t]*(-Sin[t]*
        Cos[ArcTan[a*t]] + (a*Cos[t]*(-Sin[ArcTan[a*t]]))/(1 + 
          a*t^2)), 
   y'[t] == (Sin[x[t]] - b*z[t])*
     Sin[t]*(Cos[t]*
        Cos[ArcTan[a*t]] + (a*Sin[t]*(-Sin[ArcTan[a*t]]))/(1 + 
          a*t^2)), 
   z'[t] == (Sin[z[t]] - b*y[t])*((-Cos[ArcTan[a*t]]*a)/(1 + a*t^2)), 
   x[0] == i, y[0] == u, z[0] == o}, {x, y, z}, {t, -100, 100}, {a, b,
    i, u, o}] 

I chose a,b,i,u,o as parameters, you should probably narrow that list. You can now plot results. I do not understand what you want to show exactly so I made up a plot:

Plot[Evaluate@
  Table[Tooltip[y[.1, 0.2, 0, 0, s][t] /. sol, s], {s, 0.3, 1.9, 
    0.1}], {t, 0, 10}]

The mouse will give you the parameter value. You can forego the tooltip and make a grid of plots or any other format.

Regards,

Neil