@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.

@torabi In your new boundary conditions bcs3 appear an l (lower case L, not 1). It should be 1.

Better and safer is it to include the bcs when making the variable change and indeed just make the change on dsys3:

SYS:=PDEtools:-dchange({x = L*y, Phi(x) = R3(y), psi(x) = R4(y), u(x) = R1(y), w(x) = R2(y)}, dsys3, [R1, R2, R3, R4, y]);
## But you still have the convergence problem:
res:=dsolve(SYS, numeric, output = listprocedure, initmesh = 3024, approxsoln = [R1(y) = y^2*(1-y)^2, R2(y) = sin(Pi*y), R3(y) = y*(1-y), R4(y) = y*(y-1)], abserr = 0.1e-4);

@tomleslie You are simply forgetting the "extra" r in the conversion from cartesian to polar coordinates. It should be
int( r*Heaviside(1-r^2),  r=0..infinity, theta=0..2*Pi);

@torabi I suppose that your problem has some physics or engineering origin. Thus my guess is that you should have some expectation of the way a solution could or should look. It doesn't have to be very precise, but should be of the form
Phi(x)= some concrete function of x, psi(x) = some concrete function of x, etc.

Without that, it would take me much more time than I'm willing to put into it, and success wouldn't even be sure.
That your boundary conditions for u and w seem rather complicated obviously doesn't help.

@shahid You will find it here:

http://www.maplesoft.com/applications/view.aspx?SID=1343

Download the code (a zip file).

An interesting case.
Although one should think that the following two would be handled similarly by sqrt, they aren't.
By OK I mean factored:

restart;
ifactor(149^2*157);
ifactor(151^2*157);
sqrt(149^2*157); #OK
sqrt(151^2*157);#Not OK
j:=1:
simplify(((148+j)^2*157)^(1/2));#OK
simplify(((150+j)^2*157)^(1/2));#OK
sqrt((148+j)^2*157); #OK
sqrt((150+j)^2*157); #Not OK
simplify(sqrt((150+j)^2*157)); #Not OK

@torabi From your previous questions it appears that you do know how to specify an approximate solution. See also the help page for dsolve,numeric,bvp where the optional argument approxsoln = ... is described.

I think that you mean which approximate solution would be an approiate one, i.e. will lead to convergence.So use any knowledge or expectation you have for the solution in setting up a guess (an approximate solution). I don't have any suggestion, but try yourself.