@umar khan You can use v(0) = v0 (where v0 is unassigned) and proceed as I did.
The line

res2:=dsolve(eval({sys2,sys3},res1),numeric);

  you change to

res2:=dsolve(eval({sys2,sys3, w(688240)=-2.05684},res1),numeric);

Then dsolve will with the extra condition w(688240)=-2.05684  (if possible)  find v0 as well as u and w.

You will need to use method=bvp[midrich] or method=bvp[middefer].

I had no luck with my attemps though.
 

@vv Thanks. That should do it with _MaxSols set sufficiently high.
_MaxSols:=200 should do it.

@vv Indeed, solve does a good job (here for the case m=125):

r:=solve(7/18-(1/2)*cos(15625*Pi*x)=0,x,allsolutions);
S:={op(indets(r,suffixed(_B,integer)))=b,op(indets(r,suffixed(_Z,integer)))=z};
R:=subs(S,r);
Here b = 0 or 1, z an arbitrary integer.
But he only wants roots in the interval 6/125 .. 7/125 (still for m=125).

@  I get into problems before you do.

For the record let me collect the code:
 

restart;
CreaCos := proc (C, n, m, t) local k, F; 
   F := C[1][1]+(C[1][2]-C[1][1])*t; 
   for k to n-1 do 
     F := F, C[k+1][1]+((1/2)*C[k+1][2]-(1/2)*C[k+1][1])*(1-cos(Pi*m^k*t)) 
   end do; 
   return F 
end proc;

CreaC := proc (n, c1, c2) local i, H; 
  H := NULL; 
  for i to n do 
    H := H, [c1, c2] 
  end do; 
  return [H] 
end proc;

Rdensf := proc (C, V0, n, m) local k, k0, R, raiz, V, t1, alpha; 
    alpha := evalf(sqrt(n-1)/m); 
    t1 := -10; 
    V := NULL; 
    for k to m do 
      if evalf((k-1)/m) <= evalf(V0[1]) and evalf(V0[1]) <= evalf(k/m) then 
        k0 := k; 
        break 
      end if 
    end do; 
    R := sort([RootFinding:-Analytic(CreaCos(C, n, m, x)[n]-V0[n], x, re = (k0-1)/m .. k0/m, im = 0 .. .1)]); 
    for raiz in R do 
      if evalf(Im(raiz)) <= 1/100000000 then 
        V := evalf(<CreaCos(C, n, m, raiz)>); 
        if evalf(LinearAlgebra:-Norm(V-V0)) <= alpha then 
          t1 := raiz; 
          break 
        end if 
      end if 
    end do; 
    return t1 
end proc;
##
## I run into problems when m=125:
t0 := Rdensf(CreaC(3, 0, 1), <1/18, 1/9, 1/9>, 3, 50);
t0 := Rdensf(CreaC(3, 0, 1), <1/18, 1/9, 1/9>, 3, 100);
t0 := Rdensf(CreaC(3, 0, 1), <1/18, 1/9, 1/9>, 3, 124);
t0 := Rdensf(CreaC(3, 0, 1), <1/18, 1/9, 1/9>, 3, 125);

I ran the code in Maple 2016, 2015, and 15.
In the latter case my patience ran out for the case m=100, but m=124 was OK.
All three complained at m = 125.

The call to RootFinding:-Analytic for m=125 is
RootFinding:-Analytic(7/18-(1/2)*cos(15625*Pi*x), x, re = 6/125 .. 7/125, im = 0 .. .1);
and that for me produces an error, where the corresponding command for m=124 just computes the roots.

@John Fredsted Apparently MaplePrimes was just updated.
Here is the same test today:
eq := a+b/(s+t)+c = int(f(x), x = a .. b)

So the result of using the context menu in Maple using copy special/copy as image is now that you get the editable code and the image, not quite the intention, but it may have its advantages.

Nobody answered yet. That may be for the simple reason that we cannot see the command that produced the error.
Your sentence
   " Then, for k=50, 100, 150... the instruction "
is not complete. What is the instruction?
Also I notice that you write:  "These are the procedures:"
and we only see one, CreaCos.

@umar khan In your uploaded worksheet you have not made the two changes that I suggested and which Tom Leslie seconded:

The Pi/pi problem and the exp/e problem.

That is OK, but you ought to tell us why not.

@tomleslie student_md did implicitly say that a thing like 2+3 works, so an example of something that doesn't work would be good.

@student_md Correct. Nothing is assigned to x1(t) or x2(t). Although that can be done rather easily, I would advise against it. I find it more convenient to leave x1 and x2 unassigned so they can be reused in another ode (or the same ode) in the same session without unassigning.
Here you can do:

ode2:=diff~(X,t)=A.X+B;
ic:=eval(X,t=0)=<-3,5>;
sol:=dsolve({ode2,ic}); # Giving the result a name
## The following construction may look weird the first time, but it ensures that X1 is assigned to the x1(t) result and X2 to the x2(t) result:
#
X1,X2:=op(subs(sol,[x1(t),x2(t)])); #Using capital X1, X2
# To understand it, let us split it up:
subs(sol,[x1(t),x2(t)]);
## then remove the list brackets to get a sequence of length 2
op(%);
X1,X2:=%; #Here assignment to both X1 and X2.
#Plotting
plot(X1,t=-2..5);
plot([X1,X2],t=-2..5); #Both
plot([X1,X2,t=-2..5]); # Phase plane plot: Notice the position of ].

@Kitonum I suppose that there is an error (weakness or what have you) in Maple 2016.1, which has been corrected in Maple 2016.2, the version of which I'm the lucky user.
When I reported no errors in the 4 Maple versions I was using the worksheet provided in 2D math input.

When I try your simple code in a worksheet with 1D input:

restart;
with(VectorCalculus);
int(1,[x=0..1,y=0..1]);

I get no error in Maple 2016.2, but do get errors in Maple versions 2015.2, 18, 15, and 12.

When I copy the 2D triple integral from the worksheet and plant in on an input line accepting 1D-mat input I get in Maple 2016.2:

int(Divergens, [z = rhs(S1) .. rhs(S2), y = -rhs(Bounds[1, 2]) .. rhs(Bounds[1, 2]), x = -2 .. 0]);

whereas in the other earlier versions I get

int(int(int(Divergens, z = rhs(S1) .. rhs(S2)), y = -rhs(Bounds[1, 2]) .. rhs(Bounds[1, 2])), x = -2 .. 0);

Thus the reason why it worked for me in the early  verions was that the copying of 2D math input was rendered as an int(int(int(... ))) rather than as an int(..,[...] ).

Isn't 2D-math input wonderful? (sarcasm intended).

@student_md With the inhomogeneous term and initial conditions the simplest is to use dsolve:

B:=<-t,-2*t-1>;
ode2:=diff~(X,t)=A.X+B;
ic:=eval(X,t=0)=<-3,5>;
dsolve({ode2,ic});
X=subs(%,X);
### Using the familiar formula for solving a linear first order ode, a formula that also holds for the matrix case)  can also be done:

ODE:=diff(x(t),t)+a*x(t)=b(t); #Just a reminder
dsolve({ODE,x(0)=x0});
x(t)=exp(-a*t)*x0+exp(-a*t)*Int(exp(a*s)*b(s),s=0..t); #Rewritten
In our case a = -A, and the correct order has to be respected:
E.rhs(ic)+E.Int(eval(E,t=-s).eval(B,t=s),s=0..t); #Just an exhibit!
E.rhs(ic)+E.int~(eval(E,t=-s).eval(B,t=s),s=0..t); #Now computing
combine~(expand~(%)); #Final result agreeing with the dsolve result.
## Even the first mentioned approach could be carried through. There you would first need a particular solution, which could be found by an ansatz
pol:=<a+b*t,c+d*t>;
(pol can be found to be: <(5/2)*t-13/4, -(11/2)*t+21/4>  )
Then add that to the general solution to the homogeneous equation. Finally determine the two constants c[1] and c[2] from the initial condition ic.

@Kitonum Did you actually get any errors in the original worksheet?

As I reported, I didn't get any errors in Maple 2016 (actually 2016.2), Maple 2015.2, Maple 15, and Maple 12. So I wonder why?

Your worksheet produces no errors in Maple 2015 or 2016. Not in Maple 12 or 15 either!
So which version are you using?

The last integral returns (10/3)*Pi.

Your triple integral in Maple notation (I loathe 2D input):
int(Divergens, [z = rhs(S1) .. rhs(S2), y = -rhs(Bounds[1, 2]) .. rhs(Bounds[1, 2]), x = -2 .. 0]);
Change int to Int (the inert version of int) to see the bounds evaluated:
Int(Divergens, [z = rhs(S1) .. rhs(S2), y = -rhs(Bounds[1, 2]) .. rhs(Bounds[1, 2]), x = -2 .. 0]);
## Then use value:
value(%); # result 10*Pi/3.

With your system called SYS I tried (in Maple 2016):
dsolve(SYS,numeric,method=rkf45_dae);

and got the error message:
Error, (in evala/Factors) Maple was unable to allocate enough memory to complete this computation.  Please see ?alloc

Without giving the method (or when using Maple 2015)  I got the error you reported.

I tried the context menu in Maple using copy special/copy as image
and quite right, surprisingly, I got:

eq := a+b/(s+t)+c = int(f(x), x = a .. b)

Windows 10, 64bit, and Firefox 50.0.