@Axel Vogt Your way is nice and simple!
 

rationalize(expr); 
res:=int(rationalize(expr),p); 
simplify(diff(res,p)-expr); # 0

This even works in Maple 8. In fact in Maple 8 rationalize is not necessary.
Maple 8 comes up with the following without assumptions just from res:=int(expr,p);
 

res := 1/2*(2*arctan(p)*a^2+2*arctan(p)*a^2*p^2+4*arctanh(p/(a^2-1)^(1/2))*(a^2-1)^(1/2)*p^2-((p^2+1)^2*(a^2-1))^(1/2)*ln(p^2+1)+4*arctanh(p/(a^2-1)^(1/2))*(a^2-1)^(1/2)+2*((p^2+1)^2*(a^2-1))^(1/2)*ln((p^3*(a^2-1)^(1/2)+p*(a^2-1)^(1/2)-((p^2+1)^2*(a^2-1))^(1/2)*((a-1)*(a+1))^(1/2))/(p^2+1)/(a^2-1)^(1/2))-2*arctan(p)-2*arctan(p)*p^2)/((p^2+1)^2*(a^2-1))^(1/2);
###
simplify(diff(res,p)-expr); # 0

Obviously you didn't give us the full code.

Could you do that or simply upload a worksheet with the whole shebang?

Runge-Kutta methods are used for initial value problems, but you have a boundary value problem.

@Carl Love Exploring this just one step further indicates a bug in Maple 2020.1, Windows 10, June 10 2020 Build ID 1474787.

restart;
u:=1/3*ln(_C1^2 - 2*_C1*x + x^2 + 1) + 1/2*(-2*x + 2*_C1)*arctan(-x + _C1);
simplify(u);

I get
1/3*ln(_C1^2 - 2*_C1*x + x^2 + 1) + 1/3*(-3*x + 3*_C1)*arctan(-x + _C1)

which is correct, but couldn't possibly be intended.
######################

A simpler example:
 

restart;
u:=1/3*g(a) + (x+y)*f(b);
simplify(u);

Result:
1/3*(3*x + 3*y)*f(b) + 1/3*g(a)

PS. This goes all the way back to Maple 2016.2. It is not present in Maple 2015.2.

@Carl Love I need to use simplify directly on that term:
 

restart;
interface(version); # `Standard Worksheet Interface, Maple 2020.1, Windows 10, June 10 2020 Build ID 1474787`
my_sol:=y(x) = arctan(x - _C1)*x - arctan(x - _C1)*_C1 - ln(1 + (x - _C1)^2)/2;
simplify(my_sol); # 2 is still there
map(simplify,rhs(%)); # 2 is gone

@nm dsolve/numeric/bvp comes up with the initial profile n(x) = 0, which obviously is no good since n(x) is a factor in the denominator of the expression for n''(x).

So if 'Vortex' has an idea about the shape of a solution he could try giving an approximate solution in the form
approxsoln=[n(x) = f(x)]
where f(x) is something having that shape.

@AHSAN Since you don't say what it is that needs more explanation I shall make a few comments about what I think may require some elaboration. I will restrict myself to the last and fastest version.

1. The use of unapply. In your own worksheet you have the lines:
 

lambda1 := -3*(7*k^3*sigma^3 + 32*Q*k^2*sigma^2 - 11*k^2*sigma^3 + 54*Q^2*k*sigma - 44*Q*k*sigma^2 + 11*k*sigma^3 + 36*Q^3 - 54*Q^2*sigma + 32*Q*sigma^2 - 7*sigma^3)/(20*sigma^4);
## where sigma is given in terms of x earlier.
data1 := [seq([lambda1(x), x], x = 0 .. 0.6, 0.1)];

Since lambda1 is NOT defined as a function, but used as one anyway you may wonder why that works (it does!).
To see why, try this version instead:
 

data1 := [seq([lambda1, x], x = 0 .. 0.6, 0.1)];

Notice that you get exactly the same as before because the x inside lambda1 is given the different values of x = 0..0.6 with spacing 0.1.  The first version with lambda1(x) works because any number (e.g. -0.09786536940) will also work as the constant function with that value, thus -0.09786536940(x) = -0.09786536940 for any x.

Since I'm not going to use seq in my code I will turn lambda1 into an actual function of x. That is done by unapply.
After that I can do e.g. lambda1( 0.3) and expect to get what I intended.

2.  N is chosen so that the spacing between the x-values is 10^(-6). The vector V contains all the x-values and is given the datatype float (or explicitly float[8 ] ). Otherwise it would have been datatype=anything. To see that, try:
 

V:=Vector(7,i->0.6/6*(i-1));
op(3,V);

Why do I care about the datatype? Because for speed when V is big (you will have N = 600000) I would like to use hardware float computation (evalhf) and to get the most out of that I avoid data conversion from 'anything' to 'float[8]' by setting datatype=float[8] to begin with.
3. The use of map instead of the easier elementwise operation ~ .
Try

V:=Vector(7,i->0.6/6*(i-1),datatype=float[8]);
sin~(V); # OK
map(sin,V); # OK
evalhf(sin~(V)); # error
evalhf(map(sin,V)); # OK

4. Finally < W | V > creates a matrix with the columns W and V.

Is the vector you give named gln so that gln[i], i = 1..2 are its components?

@tomleslie dsolve, events works in Maple 12, thus also in Maple 13.

Your suggestion [[ gln[1]^2+gln[2]^2, halt]]  should work in principle, but I would like to see the actual system.
Will gln[1]^2+gln[2]^2 ever be found equal to zero, problem being that it is >=0 always?

The solution to that is the obvious: Replace gln[1]^2+gln[2]^2 by gln[1]^2+gln[2]^2-epsilon, where epsilon could be e.g. 1e-7.

There is a dog chasing jogger example in the help page for stop_cond in Maple 8. stop_cond has been superseded by events.
Dog_chasing_jogger_events.mw

For your worksheet to run without error you need kernelopts(floatPi=true).
This option exists in Maple 2015, but is not documented. It doesn't exist in earlier versions. I believe the default is false in 2015.
The option is documented in Maple 2016 and later versions and the default is true.
I found out since I have kernelopts(floatPi=false) in my maple.ini file.

If a univariate polynomial has a multiple root then that root is given as many times as its multiplicity.

I suppose that if in the solving process for you equation a univariate polynomial turns up this could create the situation you describe.

Try solving with infolevel[solve]:=3.
In the lines you will see polynomials mentioned. At the end a confusing message saying:
Main: Exiting solver returning 1 solution

@ogunmiloro Which Maple release/version are you using?
More importantly, what is your response to mmcdara's question about the fact that C__f has 110 members and 'times'  has only 20? Should the members in C__f be grouped in 20 groups of varying size and so how?

I haven't as yet found a description of the meaning of the procedure _pexports.

But it appears to be short for package exports and works this way:
 

M:=module() option package; export _pexports,a,b,c;
   _pexports:=proc() 
              [op({exports(M)} minus {':-c', ':-_pexports'})]
   end proc;
   a:=proc(x) 8*x end proc;
   b:=proc(x) x^2 end proc;
   c:=proc(x) sin(x) end proc;
end module;
##########
with(M); # Notice that only [a,b] is returned
a(7);
b(s);
c(8); # unevaluated
M:-c(8); # the long form works.

It most likely goes back to the introduction of modules in Maple (Maple 6 I believe). The code above certainly works in Maple 8 and later versions.
A neat facility which makes it possible to make a useful but "nerdy" procedure available to the programmer in other contexts than within M. I'm somewhat embarrassed that I haven't noticed this before (or maybe I have, but forgot).

@dharr Yes, this extra semicolon is allowed in Maple 2019 and 2020, but not before.
Trying in those two Maple versions:
p:=proc(x); local y; y:=x; y end proc;
we notice that the parser removes the semicolon.

Here is a variant where T is defined as a procedure with option remember and returning unevaluated if t is 0 or just the global name t..

restart;

ode2 := diff(varphi(t), t, t) + omega^2*sin(varphi(t));
p0 := evalf(10/180*Pi);
te:=6:
event2 := [[diff(varphi(t), t), T(t)=t]];
##
T:=proc(t) option remember; 
  if t::identical(':-t') or t=0 then 
    'T(t)' 
  else
     subs(ld2(t),T(':-t')) 
  end if
end proc;
##
ld2 := dsolve([eval(ode2, omega = 2*Pi), varphi(0) = p0, D(varphi)(0) = 0,T(0)=0], numeric,
               discrete_variables=[T(t)::float], events = event2,abserr=1e-12,relerr=1e-9);

op(4,eval(T)); # So far the only members of the remember table
plot(T(t),t=0..te); # Creates many members
{entries(op(4,eval(T)),nolist)} minus {T(0),T(t)};

The last line returns
{0., 0.5009535940730963, 1.0019071881466945, 1.5028607822207958, 2.0038143762954004, 2.504767970370508, 3.0057215644461186, 3.506675158522232, 4.0076287525988485, 4.508582346675968, 5.009535940753591, 5.510489534831716}

This one works as is in Maple 12 too.