I believe that you're using Maple 2015 (because you've said so many times). You should always use the check boxes to put your Maple version in Question headers.

In Maple 2015, the precedence of the inert % operators is set incorrectly. This has been corrected in Maple 2020.

After you put M=0.5 in params, you get the error message

initial Newton iteration not converging

What I'm about the tell you applies only to this specific error message (and only if it starts with "initial"^[*1]). The technique is called a continuation parameter. You include a coefficient (which I'll call C here) in certain terms in the ODE system. You often need to experiment to find the right terms to use. C will vary from 0 to 1 during the initial phase of dsolve's solution process; so when C=0 it's as if the specified terms aren't there, and when C=1 ​​​​​​it's as if C wasn't there. Thus, it makes sense to include C in the complicated nonlinear terms. 

For this specific problem, I put C in every nonlinear term. Then I included option continuation= C in the dsolve command. It worked.

I also took out maxmesh=5000. Why would you impose a resource limitation (that may not even be needed) when you haven't even solved the problem yet?

[*1] Thus you need to be careful when discussing these errors here to distinguish between "Newton iteration not converging" and "initial Newton iteration not converging". Both of these occur quite often when solving boundary-layer BVPs.

The index (the part in square brackets) of a function call is accessible inside the function via the special syntax op(procname). So this works:

div5:= proc(f) local x; unapply(f(x)/5, x) end proc:
f:= x-> x^2+op(1, procname):
div5(f[5]);

The 1 is to select the first index, in case there are more than one. If there will always be exactly one, you can use simply op(procname).
 

My, what a trivial procedure, hardly worth naming, let alone writing a help page for. Anyway, here it is:

Trapz:= (X, Y)-> (X[2..]-X[..-2]).(Y[2..]+Y[..-2])/2:
trapz:= proc()
local X,Y;
    if nargs=1 then Y:= args[1]; X:= Vector([upperbound(Y)][1], i-> i)
    else X:= args[1]; Y:= args[2] 
    fi;
    if (local d:= rtable_num_dims(Y)) = 1 then Trapz(X,Y) 
    elif d=2 then Vector[row]([upperbound(Y)][2], i-> Trapz(X, Y[..,i]))
    fi
end proc:

I didn't handle the case where Y has more than 2 dimensions.

You can enter the system, numerically solve it, and plot the solution in 3D like this:

sys:= { 
    diff(x(t),t) = cos(x(t)), x(0) = 1, 
    diff(y(t),t) = sin(z(t)), y(0) = -1,
    diff(z(t),t) = z(t)*y(t), z(0) = 2
}:
sol:= dsolve(sys, numeric):
plots:-odeplot(sol, [x,y,z](t), t= 0..5);

If you want to solve it with Picard iteration, you'll need to code that yourself. Here's the plot produced by the above:

Do this:

f:= (x-1)/(x+2);
s:= convert(f, FormalPowerSeries);
c:= unapply((S-> (op(1,S),[op([2,1],S)]))(indets(s, specfunc(Sum))[]));
L:= limit(abs(c(k+1)/c(k)), k= infinity);
solve(L < 1, x); 
                   (-2,2)

Notes:

1. None of your simplify commands were needed, but neither would it be wrong for you to put them back.
2. Your line sum(k-> 1/k^2, k= 1..1000) is nonsense, and its result, 1000*f, is also nonsense. The fact that this doesn't give an error message should be considered a bug in Maple.
3. The line the turns the general term of the series into the function c is now fully automatic, no need to copy-and-paste.
4. The expression for the ratio test is abs(c(k+1)/c(k)). The extra x powers that you used are mathematical errors.
5. The final result is given in interval notation, the open interval from -2 to 2.
 

Here is how you can "automatically" add a plot to an existing plot:

AddPlot:= (p, P)-> plots:-display(P, p, _rest):
f:= (x-1)/(x+2):
p1:= plot(f, x= -3..3, discont, view= [-3..3, -5..5]):
p2:= AddPlot(
    plottools:-circle([0,0], 1, color= blue),
    p1,
    scaling= constrained
);

The following simple code is probably adequate for your needs. It could be made slightly more efficient by pre-sieving for the primes.

MyProperty:= proc(q::prime, P::name) 
local p:= ithprime(q - numtheory:-pi(q));
    if p >= q then P:= p; true else false fi
end proc
:    
MySeq:= proc(n::posint)
local k:= 0, R:= Vector(n), p, q:= 1;
    while k < n do 
        q:= nextprime(q);
        if MyProperty(q, 'p') then
            k:= k+1;
            R[k]:= [p,q];
        fi
    od;
    [seq](R)
end proc
:
S:= MySeq(9);
S := [[2, 2], [13, 11], [17, 13], [29, 17], [31, 19], [43, 23], 
  [67, 29], [71, 31], [97, 37], [107, 41]]

(p_seq, q_seq):= map(op~, [1,2], S)[];
   p_seq, q_seq := [2, 13, 17, 29, 31, 43, 67, 71, 97, 107], 
     [2, 11, 13, 17, 19, 23, 29, 31, 37, 41]

In Maple, type checking can only be done at run time. This is mostly because types can be arbitrarily complex, and checking them may involve execution of an arbitrary amount of Maple code. Of course, for the example that you show, it'd be theoretically possible before runtime; but I'm sure that you realize that your example is extremely superficial. 

The code below allows you to set the values of the parameters with sliders and have the plot update automatically. For some values of the parameters, there will be warning messages such as "cannot evaluate the solution further (left/right) of ... probably a singularity." These warnings do not prevent you from proceding with the Explore. If you want, the field arrows could be added.

restart:.
A:= a_1*c^3 - b_1:  B:= 3*a_1*c^2 - 3*b_1/c:  C:= d*c - e:
sys:= {diff(x(t), t$2) = (B*x(t)^2+C*x(t)-e)/A, x(0)=f, D(x)(0)=g}:
P:= [a_1, c, b_1, d, e, f, g]:
sol:= dsolve(sys, numeric, parameters= P):
MyPlot:= proc(P)
    sol(parameters= P);
    plots:-odeplot(
        sol,
        [x(t),diff(x(t),t)],
        t= -3..3,
        view= [-6..6, -6..6],
        thickness= 2, axes= boxed, color= red, scaling= constrained
    )
end proc
:
Explore(MyPlot(P), parameters= (P=~ -5..5));

subsindets(B, string, parse)

As Tom admits, explicitly entering the expression to be integrated is "ugly", and I find Axel's modification to be only a slight improvement. And that ugliness is made worse by the fact that the OP actually has a large number of such sets of integrals to do that I removed from the posted worksheet because the syntax issue is identical even though the integrands are not.

I like to build the expression to be integrated directly into the ODE system, like this:

restart;
eq1:=diff(f(y), y$4)+Uhs*diff(E(y),y$3)-(diff(f(y), y$2))+(diff(theta(y), y$1))= 0:
eq2:=diff(theta(y), y$2)+(diff(f(y), y$2)+1)^2+1+diff(theta(y),y$2) = 0:
E:=y->zeta*(cosh(k/2*(h1+h2-2*y)))/(cosh(k/2*(h1-h2))):
bcs:=f(h1) = -(1/2)*(Q-1-d), f(h2) = (1/2)*(Q-1-d), (D(f))(h1) = -1, (D(f))(h2) = -1,theta(h1) = 0, theta(h2) = 1:
epsilon1:=0.1:d:=1:omega:=Pi/6:h1:=-(1+epsilon1*sin(2*Pi*x)):h2:=d+epsilon2*sin(2*Pi*x+omega):F:= Q-1-d:epsilon2:=0.5:x:=1:
alpha:=Pi/6:
de:=eq1,eq2,bcs:
d1 := subs(Uhs =-2, zeta=3,k=1,{de}):
param:= {Uhs =-2, zeta=3,k=1}:
P1:= eval(diff(f(y), y$3)+Uhs*diff(E(y),y$2)-(diff(f(y), y$1)+1)+(theta(y))+sin(alpha),param):
ec:=0.5:
for Q from -3 to 3 by ec do 
    F2[Q]:= dsolve(
        d1 union {diff(J(y),y)=P1, J(h1)=0}, 
        numeric, maxmesh= 25500, continuation= lambda1
    ) 
end do:
for Q from -3 to 3 by ec do 
    P3[Q]:= eval(J(y), F2[Q](1)) - eval(J(y), F2[Q](0))  
end do;
                   P3[-3] := 2.30876985489011
                  P3[-2.5] := 1.36946865887356
                 P3[-2.0] := 0.460748411800256
                 P3[-1.5] := -0.417186862674401
                 P3[-1.0] := -1.26412891537822
                 P3[-0.5] := -2.07986512872413
                  P3[0.] := -2.86417840284371
                  P3[0.5] := -3.61684700714295
                  P3[1.0] := -4.33764445144679
                  P3[1.5] := -5.02633933429334
                  P3[2.0] := -5.68269517762163
                  P3[2.5] := -6.30647027769476
                  P3[3.0] := -6.89741755441151

You can place the cursor wherever you want (such as after the x in dx), press backspace, and change it to y. Likewise, place the cursor before the "d" to enter/edit the integrand, then start a new definite integral there. There's no need to use parentheses, but I suppose that you could if you really wanted to. In the hyperlink posted below, I have entered the expression which in standard Maple Notation would be

Int(Int(f(x,y), x= c(y)..d(y)), y= a..b)

https://learn.maplesoft.com/index.html#/?d=DQJKHPFIKPOTBIGPKTGHCFFPPTKSDPBSDGLKFNFUFIGUKNBQAFDLJOOKNHDUERMQKKCSJKPGASGFEIPQLNCUIMMRHKBFMHHGANER

I think that the problem is simply that you've made a mental arithmetic error in determing your expected result, because x^(2*k)/x^k is x^k, not x^2.

An unrelated tip: You shouldn't assume that the bound variable of the Sum is literally k (even if that's always true). Instead, extract it via op([2,1], s) or, equivalently, lhs(op(2,s)).

The code that acer wrote for you that uses operators such as #mi(...) and #mo(...) was intended to display output via Maple's GUI, not via Latex. I vaguely recall him explicitly saying that, but upon rereading, I can't find it, so maybe he didn't say it. Anyway, I'm certain that that was his intention. Those operators are part of the language HTML (I believe), and not the language Latex.

Now I wouldn't be at all surprised if Latex had some facility for interpretting HTML. If it does, it's likely possible to write a brief Maple procedure to edit Maple's Latex output in such a way that that facility was turned on.

From my personal observation of the expansion of text-formatting languages, there was a rapid divergence about 25 years ago between those primarily intended for paper documents (such as Latex and PDF) and those primarily intended for on-screen documents (such as HTML and XML). Note that I said "primarily": Almost all of them can be used on either medium, but for almost all of them there's one medium that's distinctly much easier. I'd be quite interested in reading here anyone else's opinion on this divergence.