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!

    (* 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 = NDSolve[{r''[t] == r[t]*\[Phi]'[t]^2 - 1/r[t], \[Phi]'[t] == L/
        r[t]^2, \[Phi][0] == vTangential/r0, r[0] == r0, 
       r'[0] == vRadial}, {r, \[Phi]}, {t, 0, 200}]
    (* Plotting the solution *)
    ParametricPlot[
     Evaluate[{r[t]*Cos[\[Phi][t]], r[t]*Sin[\[Phi][t]]} /. s], {t, 0, 
      2*Pi}, GridLines -> Automatic, Frame -> True]