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.

@Al86 Who knows, there might be something directly wrong in the code, although I don't think so. There is no doubt, however, that it is not very efficient.
I'd like to repeat what I said about the code in the beginning of the answer above:

"Here is another approach at solving the integral equation.
In its present version it is rather crude. I have to admit to knowing nothing about numerical solutions to integral equations. But here it is."

Anyway, why do you want to use 100 points?

Could you tell us what problem1 is?  What are the inhomogeneous terms? Better give the equations, because the homogeneous equations are not written on the form xxxxxx = 0

The roots are the same as the roots of the polynomial f1:

f1:=x^3-442*x^2+65107*x-3196058;
fsolve(f1=0,x);
    143.2030367, 148.1579896, 150.6389737

How can we help you in finding your error if you don't tell us what you did?

@roya2001 Did you try what I told you in my latest comment?

I don't understand your first question.

Use the successful indices given by
indx := indices(res, nolist); #The sequence of successful indices for results

@roya2001 Remove the definition of RES after the loop.

Add this after the loop:

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

The last sequence gives you the omega2 values for each of the successes. I got 0.95... for all with guess omega2=1.

To find another omega2 try changing the guess for omega2. In particular try omega2=2.
Before you do, insert a line clearing the table res just before the loop:
res:='res';
If you don't then previous results from another computation will be kept if the new loop has unsuccessful results where the old had successful.

With omega2=2 in the approxsoln my sequence was of length 2 and omega2= 2.33011874639388.

The next omega2 will be 4.10059501931527.

@roya2001 You did not make the change I told you about. As far as I can see you didn't change a thing.

With the change I already mentioned, it works.

@roya2001 I made he following change in your worksheet:

newsys2 := subs({omega^2 = omega2*10^13, h2(theta) = 10^(-10)*h2(theta)}, newsys);

I erased the definition of newsys22.

After that I ran dsolve as you did with newsys2 instead of newsys22.

You find omega2 by doing e.g.
res(0);

Plotting is done as usual, e.g.

plots:-odeplot(res, [theta, h1(theta)], 0 .. 1);

which should give you

@roya2001 I didn't realize till now that I was not answering farahnaz.

I suggest that you ask a separate question in another thread. Otherwise I am going to get toally confused too.
When you do, I think you ought to explain what it is you are trying to do; that is thoroughly missing.

@roya2001 The value of omega2 (scaled as above) of 0.956... corresponds to solutions for h1,h3,h4 with no zeros in the interior of the interval 0..1. For h2 there is one:




Since h2 is so small compared to the other, it makes sense to scale h2 at the same time as omega2:
newsys2 := subs({omega^2 = 10^13*omega2, h2(theta) = 10^(-10)*h2(theta)}, newsys);
I didn't do this in producing the images above, but in fact it works better: more immediate successes in the b-loop.

I suggest that you show us the failed attempt (all of it) e.g. in the form of an uploaded worksheet.