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