Basically, I have a set of differential equations that I need to solve for exactly 100 different initial conditions (given as lists for each initial condition), and then plot each solution.

Here is some sample code where I have set vrad, vtan, and deltaR (arrays of initial conditions) to an array of length two. So, given the arrays vrad, vtan, deltaR (our initial conditions) I want to be able to essentially do what this code does but for the array of solutions. Cheers!

Edit: I think I've nearly done it, I just need Table to not iterate through every tuple, but instead by index, anyone know how to do this?

    (* Scaling Quantities *)
    V = 200;
    R = 10^4;
    (* Random Quantities *)
    vrad = {0, 5};
    vtan = {0, 5};
    deltaR = {0, 5};
    (* Converting to dimensionless quantities *)
    vRadial = (V + vrad)/V;
    vTangential = (V + vtan)/V;
    r0 = (10^4 + deltaR)/R;
    L = r0*vTangential;
    (* numerical solution *)
    s = Partition[
      Flatten@Table[
        NDSolve[{r''[t] == r[t]*\[Phi]'[t]^2 - 1/r[t], \[Phi]'[t] == d/
           r[t]^2, \[Phi][0] == a, r[0] == b, 
          r'[0] == c}, {r, \[Phi]}, {t, 0, 200}], {a, vTangential/r0}, {b,
          r0}, {c, vRadial}, {d, L}], 2]
    (* Plotting the solution *)
    ParametricPlot[
     Evaluate[{r[t]*Cos[\[Phi][t]], r[t]*Sin[\[Phi][t]]} /. s], {t, 0, 
      2*Pi}, GridLines -> Automatic, Frame -> True]