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

@9009134 Sorry, I have no ideas. Maybe this problem simply doesn't have a nontrivial solution, here meant as a solution for which also F,K,Omega, besides theta, are nonzero in the interior of the interval.

You write " The objective is y2 at t =1." and yet you have z2:=sol1(0.5);
Shouldn't that be z2:=sol1(1);  or am I misunderstanding something?
This is surely irrelevant for the question you pose, but confuses me.
##
I think I got it: You are resetting the clock and starting anew, so the new time just has to run 0.5 for the total time to be 1.
I tried this piecewise version which, however, doesn't improve Hesse.
 

sysp:=diff(y1(t),t)=-piecewise(t<=0.5,u1+u1^2/2,u2+u2^2/2)*y1(t),
           diff(y2(t),t)=piecewise(t<=0.5,u1,u2)*y1(t);
res:=dsolve({sysp,y1(0)=1,y2(0)=0},numeric,parameters=[u1,u2],output=listprocedure);
Y2:=subs(res,y2(t));
obj4:=proc(u1,u2) 
   res(parameters=[u1,u2]);
   -Y2(1)
end proc;

@9009134 My doubts about the solution were expressed in my comment above:
http://www.mapleprimes.com/questions/218706-Error-in-Dsolvenumericbvp-Singularity#comment232672

entitled "Plot, approxsoln, and is it a solution?".

## It may be worth while mentioning that the result of this test:

odetest({F(eta)=0,K(eta)=0,theta(eta)=1-eta,Omega(eta)=0},SYS);

is  {0}, which means that the solution given, i.e. {F(eta)=0,K(eta)=0,theta(eta)=1-eta,Omega(eta)=0} is in fact a solution to SYS according to odetest. We should, however, remember the troublesome denominators K(eta) and Omega(eta).

@9009134 Since I'm not at all convinced about the "solution" found in the first case, I found it hard to summon interest for the second, I must admit. No reason to consider new parameters that also give problems, when you haven't successfully solved the first case.

@torabi I believe that there is no positive minimum.
Consider again the example I gave above.
I have modified it to fit your situation better:
 

restart;
ode:=diff(y(x),x,x)-omega2=0; #-1.2345 changed to 1 for convenience 

for b in [seq(10^(-n),n=0..10)] do
  res:=dsolve({ode,y(0)=0,y(1)=0,D(y)(1)=b},numeric);
  print(subs(res(0),omega2))
end do:
##Basically what goes on in the numerical solution with a parameter to be determined (here omega2) is this:
ode1:=subs(omega2=omega2(x),ode);
ode2:=diff(omega2(x),x)=0;
b:='b':
dsolve({ode1,ode2,y(0)=0,y(1)=0,D(y)(1)=b}); 

It is really just enough to look at the output from the exact version at the end:
{omega2(x) = 2*b, y(x) = b*x^2-b*x}
There is clearly no positive minimum for omega2. The nonnegative minimum exists and is 0, but it corresponds to y(x) = 0 for all x in 0..1.
I believe you have the same situation.

The syntax is
int(a(t)*b(t)+2*(diff(a(t), t))*(diff(b(t), t)), t);
But what did you expect to get out of it?