@vv Using numerical calculation of the infinite series:

restart;
S:=convert(hypergeom([2/3,2/3,2/3],[-1/3,4/3],z),FormalPowerSeries);
evalf[50](eval(S,z=99/100));

    -138.99623131713330917549791383444227380249824260753

@Damon Now the link is visible and works.

@Damon I cannot see any attachment, can you?

@farahnaz (1) My guess is that the minimal omega2 is 0.956187059957532*10^13, but it is just a guess.

(2) I wrote:
   To get the successes you can do:
   indices(res);

and that is what gives you the successful attempts in the do loop.
No assignment is made to res[b] when dsolve results in an error. By doing indices(res) you will see which b's resulted in an assignment of a table element res[b] to the table res.
Note: If res is just a name, then the assignment res[xxx]:=89; implicitly creates a table whose only index at this point is xxx and whose only entry is 89.

@farahnaz Not too surprising there are other eigenvalues. E.g. a smaller eigenvalue omega2 = 0.956187059957532*10^13.
To find it (and in general for this system) it makes sense to scale omega2 since it is so large.
The reason is that abserr is also measured against the enormous omega2. This makes solution impossible with small abserr.
So the suggestion is to replace omega2 by omega2*10^13:
Thus the new system newsys2S is:

newsys2S:=subs(omega2=10^13*omega2,newsys2);
##Now do

res:=dsolve(newsys2S union bcs union {bcs2,(D@@2)(h2)(1)=10^(-7)}, numeric,approxsoln = [omega2 = 1, h1(theta) = -6*theta*(1-theta), h2(theta) = 1e-5*theta*(1-theta)(1/2-theta), h3(theta) = theta*(1-theta), h4(theta) = -theta*(1-theta)], abserr = 1e-4,initmesh=256);

##You find the new eigenvalue here
res(0);

#Of course you have to remember to multiply the value by the scaling factor 10^13.

##The loop now works on the system with scaled omega2. Make the scaling change in the code before the loop. Use Digits=15. Run the loop using abserr=1e-3. In approxsoln change the omega2 guess to omega2=1 (to start with).
To get the successes you can do:
indices(res);
I found 3:
[(D(h2))(1)], [(D(h2))(0)], [((D@@2)(h1))(0)], [((D@@3)(h1))(0)]
To find the corresponding omega2's do
res[(D(h2))(1)](0); #Gives the eigenvalue found earlier 2.07982729632676 (*10^13)
res[(D(h2))(0)](0); #The same
res[((D@@2)(h1))(0)]; #omega2=6.27005104177759 (*10^13)
res[((D@@3)(h1))(0)]; #omega2= 12.5461315012613 (*10^13)
#By changing approxsoln, inparticular the value of omega you will find more eigenvalues.

@farahnaz I think that I already gave the answer to your question 2 in my comment titled "omega".

I don't quite understand question 1. But as I have emphasized in the same comment above: boundary value problems are difficult. I cannot give you a ready made solution which will work for all choices of your parameters.

@farahnaz Do 2 things:

1. Change Digits:=20: to Digits:=15:
  This makes the whole thing go much faster.
2. Add the lines I gave you above at the very bottom. These were the lines:

sol:=dsolve(newsys2 union bcs union {bcs2,(D@@2)(h2)(1)=1}, numeric,approxsoln = [omega2 = 2.08*10^13, h1(theta) = 5*10^7*theta*(1-theta)*(1/2-theta), h2(theta) = 0.05, h3(theta) = 7*10^6*theta*(1-theta)*(1/2-theta), h4(theta) = 5*10^6*theta*(1-theta)*(1/2-theta)], abserr = 0.1,initmesh=1024);

plots:-odeplot(sol,[seq([theta,cat(h,i)(theta)],i=1..4)],thickness=3);

plots:-odeplot(sol,[theta,h2(theta)]);

### Finally, run the whole thing from the restart at the top. Ignore the error messages comming from the loop.
### Once you see the bottom lines working, i.e. when you see that pictures are actually produced, then start cleaning up the mess and that could involve removing the entire loop.
You will want to keep the definition of newsys2, bcs, and bcs2 since they are used at the bottom.

@acer Responding to

But I only know of one way to make such a policy be implemented by the administrators of this site: we members (who currently have the ability to do so) would all have to stop deleting the spam.

I'll give it a try.

@farahnaz The value of omega2 (remember that we replaced omega^2 with omega2) is found by dsolve. I didn't make it up.
What I did do was to give dsolve a start value for omega2. That start value omega2 = 2.08*10^13 was found from using my homemade procedure, which requires an approximate solution.
As an approximate solution in my procedure I actually used the one we previously used:

[omega2 = 1, h1(theta) = theta*(1-theta), h2(theta) = sin(Pi*theta), h3(theta) = theta*(1-theta), h4(theta) = theta^2*(1-theta)^2]

i.e. we were trying to find a solution with zeros only at the end points.

But my procedure found a solution which more resembles the one I subsequently used on dsolve:

[omega2 = 2.08*10^13, h1(theta) = 5*10^7*theta*(1-theta)*(1/2-theta), h2(theta) = 0.05, h3(theta) = 7*10^6*theta*(1-theta)*(1/2-theta), h4(theta) = 5*10^6*theta*(1-theta)*(1/2-theta)]

i.e. a solution with zero(s) also in the interior of the interval.

When I used my own procedure I didn't use a loop nor did I use a try... end try construction. I started with the b found in our previously successful attempt for another (but somewhat similar) system, b=(D@@3)(h2)(1).
My procedure has to solve the total set of boundary conditions for (in this case) 14 unknowns. It couldn't do that when using (D@@3)(h2)(1) so I quite simply tried using (D@@3)(h2)(1) = 1 as the extra condition and had success.
End of story.

My advice is to play with these things. Don't start out with a loop etc. Try one extra condition at a time. Use all the options you know of, try changing them etc.
Boundary value problems are definitely more difficult to deal with than initial value problems.

@Christopher2222 I delete quite a few every day. It is irritating though.

@farahnaz Since you have a different system (again) I think you ought to start a separate thread for that one. That way other people will be more likely to participate. When a new problem is introduced in comments people not participating in the discussion from the start will tend to ignore it.

@Carl Love Correcting the obviously erroneous abserr in the worksheet to 1e-4 in the latter part cures the problem while keeping middefer. 
In the first part I used abserr = 3e-4 and maxmesh = 1024 instead.

You are not using an abserr of 10^(-8), which I believe you think.
If I remember correctly, I have pointed this out earlier.
The abserr you are actually using is abserr = .466507380209733.
What you are writing when you write 1.1^(-8) is in fact (1.1)^(-8), which is .466507380209733.

I don't see any problem with your worksheet, but why do you use output=Array(.....) ?
Why not leave that out entirely?
You will then have to be more careful with the range argument in odeplot. With output=Array it is ignored.

To get the full picture (without output=Array) use the range 0..10  (N=10).
To get a close-up at t=0 use e.g. 0..0.1.

@farahnaz I don't see any way that you can do this in Maple (short of coding it all yourself).

1. The boundary conditions are given on all 4 sides of a rectangle: (0,z), (x,h/2), (x,-h/2), (L,z).
Thus the requirement for pdsolve that the system be an initial and boundary value problem is not met. pdsolve uses a time based solver. Thus you can give initial conditions for time = t0 (time will be either x or z in your case) and boundary conditions on two sides perpendicular to the time axis.
2. The problem reported by the error message in your latest worksheet must be corrected: On the line t=t0 (whatever t is, x or z) you can only have terms with that t fixed (at t=t0) not the space variable. Similarly on a line with fixed space variable you can only have conditions with that fixed.
3. I see no solution to the problem I mentioned in my previous comment and exemplified.

Rest assured that even if you correct 2 out of the 3 problems above the third will prevent pdsolve from working.
pdsolve just reports the error it happens to find first and then it stops.