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

The immediate problems that you are facing as pointed out by Tom Leslie are the proper installation of Douglas Meade's Shoot library from MapleSoft Application Center.
Once that is done correctly and once you know exactly, where you installed it, you also need to replace ODE by ODE1, where
ODE1:=solve(ODE,indets(ODE,specfunc(diff)));
the point being that ODE1 has the derivatives all isolated on the left.
In ODE you have this equation (the last):
diff(fppp(eta), eta) = N1*(3.*fpp(eta)+(eta-2.*f(eta))*fppp(eta)-2.*N2*N2*m(eta)*(diff(mp(eta), eta)))

In that equation a differentiated function appears on the right.

I had success with this call:
shoot(ODE1,IC,BC,FNS,[beta1=.1,alpha1=.1,beta2=.1,beta3=1,alpha2=1]);

@John SREH If you take a look at eq, then you see three terms starting with fourier.

The last one fourier(y(t),t,w) we must consider the unknown function (agreed?), but the two other fourier terms are not directly related to Y(w)=fourier(y(t),t,w).
Do you have an ode in Y(w)? That is not at all obvious, even if it were true.
There is no way you are going to get dsolve to succeed with that equation.
####
You will have success with equ, though:

sol := dsolve({csi, equ}, numeric);
plots:-odeplot(sol,[t,y(t)],-1.4..2.0367);

It is a nice image, but what are we supposed to do with it? Give us code, not images.