This must be a solve problem. dsolve has to use solve.
Using the option implicit=true (or just 'implicit') produces the correct implicit result.
restart;
res:=dsolve({ diff(y(t),t) = y(t)^2 - y(t)^3, y(0) = 1/10 },implicit);
eval(res,{t=0,y(t)=1/10});
odetest(res,{ diff(y(t),t) = y(t)^2 - y(t)^3, y(0) = 1/10 });
plots:-implicitplot(eval(res,y(t)=y),t=-10..20,y=0..1,gridrefine=1);
##Trying to solve for y:
res2:=solve(res,y(t));
evalf(eval(res2,t=0));#Wrong
dsolve({ diff(y(t),t) = y(t)^2 - y(t)^3, y(0) = 1/10 });
evalf(eval(%,t=0)); #Same thing

#I didn't experiment enough to ever get the correct result you mention. I only got the RootOf-expression

####Continuing with the implicit result to get the correct one:
combine(res) assuming y(t)<1,y(t)>0;
evalc(%) assuming y(t)<1,y(t)>0;
res3:=solve(%,y(t));
plot(res3,t=-10..20);

@Adri van der Meer He is referring to the loop beginning with for b in extra_bc.

That loop doesn't provide a successful index, i.e. a successful extra boundary condition as it is. As you say, he already has a success before the loop.

The loop will catch that success too, if b=1 is replaces with b=10^(-7) and initmesh is set to 256.

Note. If the scaling of f2 by 10^(-7) is removed and initmesh is set to 512 then he gets 7 successes. The corresponding scaled omega2 values are 6.05330544938177 for 5 of them and 3.16869987517702 for the remaining 2.
 

@Axel Vogt Thanks for providing Maple code for Brent's method. An alternative to fsolve is needed occasionally.

@tomleslie
Just a comment.
When I try your nice worksheet with method=branchandbound in NLPSolve I get the A value
199.961306129991433460308124778
with a minimum of approximately 8e-30. It takes longer though.

When can we expect to see something done about the spam?

Is it actually being worked on by somebody?

From looking at the image I don't see any syntax errors. But I'm surprised that pdsolve returns unevaluated as if in that session of yours pdsolve was not known to the system.
(The error messages later are due to the fact that pdsolve returned unevaluated).
The spatial range given doesn't agree with the implied range r=0.00035..0.00040, and that should in fact result in an error.
Take this example from the help page:
?pdsolve,numeric
PDE := diff(u(x,t),t)=1/10*diff(u(x,t),x,x);
IBC := {u(x,0)=1, u(0,t)=0, D[1](u)(1,t)=0};
pds := pdsolve(PDE,IBC,numeric,range=-1..1); #I added a larger range than the implied 0..1

Error, (in pdsolve/numeric/par_hyp) values in range argument must match input boundary conditions

It is much easier to spot problems in an actual worksheet rather than in an image.

@Carl Love Maybe I misunderstand you, but I think it is the other way around. To find out I did the following (very unwise).

Note: The library archive BVPtools was created by me.

unprotect(solve);
solve:=proc() print("The answer is: "); 89 end proc;
savelibname:="E:/BVPtools/BVPtools.mla";
savelib(solve);

##### In a fresh worksheet I did:
libname:="E:/BVPtools",libname;

solve(x^2=1,x);
                       "The answer is: "
                               89
## I shall remove that definition of solve right away, although it returns rather consistent results.

@Carl Love Why not prepend? 
Actually, I have always used this version:

libname:= "C:/Maple_packages", libname;

So what could go wrong?

@roya2001 I apologize for overlooking the code. Anyway, remove everything after that. You see 2 successful indices.

One of your equations is rather big and contains some terms that (maybe) extremely small. The expression
tanh(3.711056549*10^(-13)/T_Bulk^2)
has a value that is approximately 0.0000015 at the GUESS value.
I tried applying fnormal to the system. With Digits=10 the reduction in size was dramatic.

systotal:=subs(NBT=NBBT,{ueq2,Teq,eq3,eq4,eq5,eq6,eq7,U(0)=0,UF(0)=0,UT(0)=0,UFT(0)=0,((-cbf*rhobf+cp*rhop)*UF(1)+ rhobf*cbf*U(1))/10000=p2,(((-cbf*rhobf+cp*rhop)*UFT(1)+ rhobf*cbf*UT(1)))/(p2*10000)=T_Bulk,UF(1)=Phiavg*U(1),D(u)(0)=0,u(1)=-lambda*D(u)(1),D(T)(0)=0,phi(0)=phi0,D(T)(1)=1}):
indets(systotal,specfunc(tanh));
eval(tanh(3.711056549*10^(-13)/T_Bulk^2),GUESS);
##Splitting:
sys,bcs:=selectremove(has,systotal,diff):
sys:=fnormal(sys);
#No immediate result here:
res:=dsolve(sys union bcs,numeric,approxsoln=GUESS,method=bvp[midrich]);

However, it should be easier to analyse the new system for problems.
One thing is noteworthy, T_Bulk doesn't appear in (the new) sys. It does appear, however, in bcs.
But as it appears there it will really just be determined from the rest of the unknowns. Thus the condition involving T_Bulk can be left out.

@Al86 fsolve has to deal with 91 nonlinear equations in as many unknowns. fsolve can be very slow in finding solutions even for much fewer equations. One reason is that it is very demanding in making sure that it really has a solution. That in itself is not a bad thing.
In this case, as opposed to the discrete method above, the equations are all big and cannot be solved from the top and down.
The total 'length' of eqs when n=91 is
length(eqs);
   6526973

The length is a measure of the size of the problem:
?length

### Finally, it is not always (never?) an advantage to have many terms in such an expansion even when one is willing to deal with the added complexity. Take the following simple example:

restart;
eq:=y(t)=sin(Pi*t);
T:=[seq(0..1,0.1)];
Y:=t->add(c[i]*t^i,i=0..10);
eqY:=eval(eq,y=Y);
eqs:={seq(evalf(eval(eqY,t=t0)),t0=T)};
fsolve(eqs);
sol:=eval(Y(t),%);
plot(sol-sin(Pi*t),t=0..1);
#This finds the error:
numapprox:-infnorm(sol-sin(Pi*t),t=0..1);
4.418620553*10^(-9)
#######
T:=[seq(0..1,0.01)];
Y:=t->add(c[i]*t^i,i=0..100);
eqY:=eval(eq,y=Y);
eqs:={seq(evalf(eval(eqY,t=t0)),t0=T)}:
length(eqs);
fsolve(eqs):
sol:=eval(Y(t),%);
plot(sol-sin(Pi*t),t=0..1);
numapprox:-infnorm(sol-sin(Pi*t),t=0..1);
2.296156998*10^(-8)

You see that actually the error is larger with n=100. OK you could increase Digits, so maybe?

The following version of the "simple minded approach" uses the trapezoidal rule instead of a simple Riemann sum.

It has also been changed somewhat, but the idea remains the same. For comparison both integration methods are shown in use. The trazoidal rule is a second order method in d whereas the Riemann sum method is first order.

Since for this integral equation the discrete equations can be solved from the top I have done that. Although replaces one call to fsolve by N calls, the equations are considerably simpler and the whole think is much faster.

restart;
eq2:=y(t)=1-h*Int(JacobiTheta3(0,exp(-Pi^2*(t-s)))*y(s)^4,s=0..t); #Better

#The discretized integrand multiplied by the length of the (equidistant) partition of 0..1 into N intervals.
#This is used as is by the Riemann sum method:
tm:=d*JacobiTheta3(0,exp(-Pi^2*(t[i]-t[j])))*y[j]^4;
#The trapezoidal term is instead:
tmTrap:=(tm+eval(tm,j=j-1))/2;
eqTrap:=y[i]=1-h*Sum(tmTrap,j=1..i-1); #Trapezoidal
eqR:=y[i]=1-h*Sum(tm,j=1..i-1); #Riemann
h:=0.65e-4;
N:=200:
d:=1./N;
T:=Vector(N+1,i->(i-1)*d,datatype=float): #The partition
S:={seq(t[i]=T[i+1],i=0..N)}: #Substitution equations
### Simple Riemann sum
t0:=time():
y:='y':
y[0]:=1:
for i to N do
   eq1:=eval(eqR,{Sum=add} union S);
   assign(fsolve(eq1,{y[i]=y[i-1]}))
end do:
time()-t0;
YR:=Vector(N,i->y[i],datatype=float):
plot(T[1..N],YR); pR:=%:
### The trapezoidal rule
t0:=time():
y:='y':
y[0]:=1:
for i to N do
   eq1:=eval(eqTrap,{Sum=add} union S);
   #eq1:=eval(value(eqTrap),S);
   assign(fsolve(eq1,{y[i]=y[i-1]}))
end do:
time()-t0;
YTrap:=Vector(N,i->y[i],datatype=float):
plot(T[1..N],YTrap,color=blue); pT:=%:
plots:-display(pR,pT);
##There is a visible difference between the results even when using N=200.

@roya2001 You didn't do what I wrote above. I told you to use the results from a call to indices. I even gave the code. Yet that code is nowhere to be seen in your latest version.
The code that has to be appended after the loop is

indx := indices(res, nolist); #The sequence of successful indices for results
nops([indx]); #I got 6 when the guess is omega2=1 (and 2 when omega2=2)
res[indx[1]]; #Example of how you fish out the list of procedures
seq(subs(res[indx[i]](0),omega2(0)),i=1..nops([indx]));

Forget about the rest; it is a remnant of earlier attempts.

Since this is obviously a student project, I think that by now I have helped you more than your instructor or professor would appreciate.

You ought to upload a worksheet. What we see are images. Use the fat green arrow in the editor window.

@darkspin1 Could you please state the problem mathematically and not by your apparently failed attempt to do it in Maple?
By that I mean: What are the inhomogeneous equations and what are the initial or boundary conditions, or what have you?

Since you are trying to piece together solutions to the left and right of xi it might be worth your while to state the problem clearly. In fact maybe piecewise could be used.