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

Taking a look at the procedures BesselJ and `evalf/BesselJ` it appears clear that the problem is in `evalf/BesselJ`.

To see the procedures do

showstat(BesselJ);
showstat(`evalf/BesselJ`);

In the first, BesselJ, you see from the first 5 lines:

showstat(BesselJ,1..5);

that what is returned from BesselJ(-3.0,0) or BesselJ(-3,0.) is in fact `evalf/BesselJ`(-3.0,0);

Now look at the lines 7..13 in `evalf/BesselJ`:
showstat(`evalf/BesselJ`,7..13);
There is a problematic statement that says if x=0 then
res := `if`(v = 0,1.+`if`(type([v, x],('list')('float')),0,0.*I),x^v);

In our case v:=-3.0 and x:=0. , so since v<>0, res evaluates to 0.^(-3.0) , which evaluates to Float(infinity).

Note: It should be pointed out that the result of evalb(0=0.); is true

@Kitonum The extraneous lines disappear (of course) when you add the option 'style=point' to the plot P.

We can conclude that here plot connects points which we don't want connected.

If you change the parametrization so it makes a continuous curve, e.g. by going counter clockwise like the following then everything should be fine.

Sq:= t-> piecewise(t < 1, [t,0], t < 2, [1,t-1], t < 3, [3-t,1], [0,4-t]):

@rlopez One looks in vain in the help for type for the type 'known'.

It is known (no pun intended) to the system from the start, however.
`type/known`evaluates to And(function, `PDEtools/known_function`).

restart;
`type/known`;
showstat(`type/PDEtools/known_function`);
## diff is considered known:
type(diff(y(x),x),known);

@catrapato In one procedure of mine I used the following method here exemplified with your system:

restart;
odex[1] := cos(x)*(diff(y[1](x), x, x))+(K^2)*sinh((1/10)*x+1)*y[2](x)+y[1](x)+exp(-x)-1:
odex[2] := (diff(y[2](x), x, x))/(1+x)+(K^2)*exp((1/10)*x+1)*y[1](x)/(1+x)-y[2](x)+erf(-x)-1:

vars:=remove(has, indets(indets({odex[1],odex[2]},specfunc(diff)),function), diff);
x:=op(op~(1,vars)); #Finding the independent variable

## The method proposed by Robert Lopez works at least as well:
Q:=indets(indets({odex[1],odex[2]},'specfunc'('anything','diff')),'And'('function','Not'('known')));

Presumably there is one more ode, one involving y[2] ?

In any case the given ode is taken out of context if y[1] and y[2] are both supposed to be unknown.

@Markiyan Hirnyk Presumably there is one more ode, one involving y[2].

In any case the given ode is taken out of context if y[1] and y[2] are both supposed to be unknown.

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