Combining "Do" and NDSolve

I have following code:

Two differential equations,each one of them with varying initial conditions

Clear["Global*"]
sol1[l0_] :=  NDSolve[{D[r[x]/(1 + x)^4, x] == 0, r[l0] == 2}, r, {x, 0,1000}, AccuracyGoal -> 15];

sol2[l0_] := NDSolve[{D[n[x]/(1 + x)^4, x] == (r[x] /. First@sol1[l0]), n[l0] == 2}, n, {x, 0, 1000}, AccuracyGoal -> 15];

Next I try to plot the solution of the sum of solutions of differential equations as function of $$l0$$, but this doesn't work:

len = 3; sum := 0;
sump[l0_, i_] := n[i] /. sol2[l0];
Do[sum += sump[l0, i] , {i, len}];
Plot[sum, {l0, 1, 2}, PlotRange -> All, PlotStyle -> Red]

Some alternative solution: If I choose "l0" in the first differential equation it works

Clear["Global*"]
sol1 := NDSolve[{D[r[x]/(1 + x)^4, x] == 0, r == 2}, r, {x, 0, 1000}, AccuracyGoal -> 15];
sol3[l0_] := NDSolve[{D[n[x]/(1 + x)^4, x] == (r[x] /. First@sol1), n[l0] == 2},n, {x, 0, 1000}, AccuracyGoal -> 15];

len = 3;
sump[l0_, i_] := n[i] /. sol3[l0];
sum := 0;
Do[sum += sump[l0, i] , {i, len}];
Plot[sum, {l0, 1, 2}, PlotRange -> All, PlotStyle -> Red]
• Have you tried evaluating sump[l0, 1], for instance? Do you notice the problem? If you do not provide a numerical value for l0, NDSolve cannot proceed symbolically and throws (many, many) errors. Jun 9 '20 at 20:48
• You are right, sorry didn't I see it before. There is some way to plot the sum of some "n/.sol[l0]+n2]/.sol[l0]+...n/.sol[l0]" for different values of l0? Jun 9 '20 at 20:56

Note that all your equations are solvable analytically:

sol1 = DSolveValue[{D[r[x]/(1 + x)^4, x] == 0, r == 2}, r[x], x];
sol3 = DSolveValue[{D[n[x]/(1 + x)^4, x] == sol1, n[l0] == 2}, n[x], x];
sum = Total[sol3 /. {x -> {1, 2, 3}}];
Plot[sum, {l0, 1, 2}, PlotRange -> All, PlotStyle -> Red] • I always try to put the simplest equations that I find in order to don't post my homework. Jun 9 '20 at 21:38
• @Cruz Yes, I figure that this might be a toy example. Look at ParametricNDSolve though: that allows you to work with a numerical solution of a differential equation that contains a parameter. You can then substitute numerical values of the parameter and carry out the calculations on multiple values when you need to. Jun 9 '20 at 21:41
• Thank you for your suggestions. Jun 9 '20 at 21:42

The suggestion of @MarcoB was really helpful.

This is my propose of solution:

Clear["Global`*"]
sol1[l0_] :=  NDSolve[{D[r[x]/(1 + x)^4, x] == 0, r[l0] == 2}, r, {x, 0,1000}, AccuracyGoal -> 15];

sol2[l0_] := NDSolve[{D[n[x]/(1 + x)^4, x] == (r[x] /. First@sol1[l0]), n[l0] == 2}, n, {x, 0, 1000}, AccuracyGoal -> 15];

len = 3;

Plot[Sum[n[i] /. sol2[l0], {i, 1, len}], {l0, 1, 2}]

Sorry for posting this question, it was my confusion.