If you haven't already, you should look at the help page for DifferentialGeometry:-JetCalculus:-TotalDiff. I think that will give you a start.

@danortega I need you to tell me any lower or upper bounds that you know of for *all* of the variables and constants:

phi__1, phi__2: (nonnegative?)
gamma: (in interval 0..1?)
pi__w, pi__a, pi__b: (in interval 0..1?)
w: (greater than 1?)
beta: (nonnegative?)

Are any of the a, b, c, d necessarily greater than 1?

In this case, you should post the original mathematical problem that you're trying to solve, rather than something (such as the code in your Question) that is at least 2 steps removed (i.e., levels of abstraction) from that original problem:

• Level 1 abstraction: You're trying to apply some symbolic solution technique (i.e., algorithm) to the problem. I suspect that the problem is a differential equation, and you're trying to find a series solution. But what you've posted doesn't contain any differential equation, nor does it contain any series.
• Level 2 abstraction: You've tried to encode Level 1 into Maple. Tom Leslie has interpreted your Question's code into working Maple code, and I would've interpreted it exactly as he did, but I doubt that it's relevant to your actual problem because it has a trivial and unique solution, given below.

The trivial solution: 

U[k](x) = 3^k/k! * exp(x)  #for all nonnegative integers k
             

Sum(U[k](x), k= 0..infinity) = exp(3+x)
              

@Carl Love 

I am surprised that dsolve does return a solution for cases (such as yours) where the delay itself depends on one of the functions being solved for; however, for the cases that I've tried, the solutions don't seem correct (using Maple 2019 at the moment). For example,

ode:= diff(y(x),x) = y(1-y(x)):
sol:= dsolve({ode, y(0)=1}, numeric, delaymax= 99, delaypts=10^5):
plots:-odeplot(sol);
 

A close inspection of that plot's data matrix shows that the plotted function is y = 1+x down to the last decimal place. Direct substitution of that into ode shows that that isn't correct.

So, is this a bug, or am I not using the delay DE solver correctly?

Here are some minor steps, mostly to simplify the presentation of the problem. I don't think that they'll have a major impact towards obtaining a symbolic solution (but every little bit helps): 

1. Please get rid of the subscript t. It's just clutter that makes the problem harder to read.
2. For the same reason, replace 1-gamma with some other variable, say g.
3. Substitute 1 - phi__1 - phi__2 for phi__3 in the objective.
4. Step 3 eliminates all decision variables from the fourth term, so that term can be removed entirely; it can't change the maximizing values of the decision variables.
5. Replace a__1 + 1 by A__1, a__2 + 1 by A__2, b__1 - 1 by B__1, b__2 - 1 by B__2, c__1 + 1 by C__1, c__2 + 1 by C__2, d__1 - 1 by D__1, d__2 - 1 by D__2, and w+1 by W.

With these changes, the objective can be entered like this:

phi:= <phi1,phi2>: A:= <A1|A2>: B:= <B1|B2>: C:= <C1|C2>: E:= <D1|D2>:

One way:
Ob:= (pi__w*((W-<2|2>.phi)^g+beta*(pi__a*(A.phi-1)^g+(1-pi__a)*(C.phi-1)^g))
    +(1-pi__w)*beta*(pi__b*(B.phi+1)^g+(1-pi__b)*(E.phi+1)^g))/g;

Or, using more-Maple-specific operators:
Ob:= (<map(`.`, [-<2|2>, A, C, B, E], phi) + [W,-1,-1,1,1]>^%T)^~g
    . <pi__w*~(1,beta*~(pi__a,1-pi__a)),(1-pi__w)*beta*~(pi__b,1-pi__b)>/g; 

Now here's what may allow some major progress towards a symbolic solution: I suspect that there are some bounds for your variables and constants that you haven't stated. Which are necessarily nonnegative? Which are necessarily in the interval 0..1?

@rcorless In addition to the errors that you've correctly pointed out, the code also has both k and k[1] used as bound^[*1] variables of summation. So k is being used four ways.

[*1] Bound here is the adjective, the opposite of free and the past participle of to bind. It's not meant in its noun sense related to the boundaries or limits of summation. 

@tomleslie Good Answer; voted up. It can be made a bit simpler because Unit is pre-defined as a top-level command. It does pass directly to Units:-Unit, but you don't need to explicitly invoke the Units package to use it.

@pik1432 The code that I gave must be done before you use with(DEtools). (That's why I said "Then proceed with the code that you already had"---though I do realize that I should've stated that more clearly; so, sorry about that.)

Since your curve is not a function in the usual mathematical sense (exactly one y-value for every x in its domain), it's not clear what you mean by "area under this function".

@pik1432 

DEtools[diffop2de] was not written to work on equations. This is a simple issue of syntax rather than a mathematical reason why it doesn't work. So, it's very easy to overload it so that it will work:

unprotect(DEtools):
DEtools[diffop2de]:= L-> 
    `if`(L::algebraic, `DEtools/diffop2de`, curry(map, procname))(args):
protect(DEtools):

Then just proceed with the code that you already had.

What is the origin of this problem? It puzzles me that someone took the trouble to come up with the interesting coefficients 1/3, 1/9, 2/3, 4/9 even though they play no role in the solution.

And why do you want to use RK2 specifically? This problem too trivial to highlight any differences between methods. Most numerical ODE solvers (including all the RKs) start by algebraically solving the system for the highest-order derivatives. In this case, that gives

diff(A[0](t), t) = 2*A[2](t),
diff(A[1](t), t) = 0, #for all t
diff(A[2](t), t) = 0 #for all t

@Rouben Rostamian  The three minima that she^[*1] is referring to are the critical points of the restriction of the objective function to a steepest-descent line, not critical points of the objective itself. Using the initial point (x,y) = (-1,1), this restricted objective is P:= t-> 16*(2*t - 1)^2*(1600*t^2 + 1). This has three critical points: two minima and one maximum. So, yes, the statement that there are three minima is indeed incorrect, but I don't think that it's incorrect in the way that you were thinking. There is still a need to pick one of them.

[*1] I guess that the OP is a woman based on the name Zeineb, which is a female name of Arabic origin.

@Rouben Rostamian The objective function here is the Rosenbrock function (see Wikipedia link), which is a standard test case for numerical minimization techniques. It's intentionally designed to have a unique local (and also global, but that's not as relevant) minimum that's trivial to compute symbolically "by hand" yet difficult numerically, especially for gradient-based techniques. Thus it's clear to me that the OP is attempting to code and test a numerical technique and is not a student taking a first course in multivariable calculus (where that "by hand" method is taught).

@Axel Vogt Sorry, I wasn't trying! to correct you. My interpretation was wrong also. I now see that the OP has the confusing habit of putting the expressions to be input on the right sides of equations and the desired output on the left. The OP also consistently mentions the desired output before mentioning the input in every sentence where both are mentioned. Both of those habits are the reverse of the usual writing style. So, I asked "Do you have it reversed?" But it turns out that you had the direction of conversion that the OP intended. It just needs to be converted from regular BesselJ to SphericalBesselJ, which seems pretty easy, although we'd need to define SBesselJ as a new symbolic mathematical function.

@Carl Love I can't think of any good reason to not split that 3rd subsequence in two: Replace f[3] with

f[3]:= n-> 4*n-1:
f[4]:= n-> 4*n:

That should make the analysis easier. The fact that these two appear so close to each other on the graphs is not a good reason to keep them together when it's obvious that 4-cycles are key to understanding the overall sequence.