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

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.