Why not use LinearAlgebra:-Eigenvector?

Choosing y ' (0) = 0 as suggested by tomleslie works fine numerically on the interval 0..1:

restart;
edo := -b*y(x)*(diff(y(x), x))^2*x-a*y(x)*x+(diff(y(x), x, x))*x+2*(diff(y(x), x)) = 0
edo1 := subs(a = 1, b = 1, edo);
res:=dsolve({edo1,y(0)=1,D(y)(0)=0},numeric,method=bvp[midrich],range=0..1);
plots:-odeplot(res,[x,y(x)],0..1);

## To push it further, you have to work harder:

res:=dsolve({edo1,y(0)=1,D(y)(0)=0},numeric,method=bvp[midrich],range=0..2,approxsoln=[y(x)=1-0.3*x^2],initmesh=256);
plots:-odeplot(res,[x,y(x)],0..2);

The solution to the singular initial value problem {edo1, y(0)=1, D(y)(0)=0} probably tends to infinity at some finite x0, i.e. y(x)-> infinity for x-> x0 from the left. x0 seems to be greater than 2 and probably less than 3.

Note. There is some evidence for an analytic solution at zero having radius of convergence less than 2.141580437. Powers all even. Using the terms up to degree 50 of that solution (sol) as approxsoln in dsolve you can get success with
res:=dsolve({edo1,y(0)=1,D(y)(0)=0},numeric,method=bvp[midrich],range=0..2.09,approxsoln=[sol],initmesh=1024);
plots:-odeplot(res);

@ecterrab We agree that a syntax like f(a)(b) certainly is useful as in e.g. D(f)(b).
I was deliberately vague when calling f(t)(t) a weird construction.

But notice that it was the appearance of several "weird" constructions like that that led me to find Torleif's mistake, which he admits to having made, when making last minute changes to a worksheet.

So far nobody asked you why you need this function. What do you intend to use it for? The answer to that may determine which method is preferable.

@ecterrab Sure, but can you give a non-weird example of f(t)(t)? Point being that there are 2 t's.

D(sin)(sin) would not make much sense, would it?

@torabi What sense does it have to remove nonlinear terms?

@torabi I had a nagging suspicion that your change to b=0 was not a good idea: If the system has the solution s=0, g=0, no what the values of omega11 and omega22 are,  then using b=0 would most likely lead to that solution.

I didn't get down to checking. But now by doing

plots:-odeplot(res[indx[i]],[[x,s(x)],[x,g(x)]],0..1,thickness=3);

for i = 1, 2, 3, ..., 11 (one at a time), you will see that 1,6,7,9,10, and maybe 11, are probably (s,g)=(0,0) with numerical roundoff.

So I suggest changing b=~0 to b=~1.

@torabi You already have them in your latest worksheet where you used the code



If you want a tag put on them, you can do this instead::

seq(subs(res[i](0),[i,omega11(0),omega22(0)]),i=indx);

Your problem is nonlinear, so I don't see any possibility for an analytical solution to the problem.

@torleian You say "I'm sorry, but I don't quite understand your earlier objection".
All my objections have been about the point that something like f(t)(t) is strange, although not syntactically wrong in Maple.

@torleian Let me state first that I'm not a user of the Physics package.
I find the objectionable term in g0. I still find it meaningless, to put it bluntly.

So to me it is no wonder that the system must object at some point. Your procedure diff_dI has to differentiate g0 w.r.t. theta(t) (and also w.r.t. diff(theta(t),t)) using the Physics version of diff.

I tried the two loops below. You get the error messages we have seen, but are the results (where there is no error) as you expected?
I find all of these expressions in IDTS below strange:
(cos(y(t)+phi(theta(t))))(t), (diff(theta(t), t))(t), (diff(theta(t), t, t))(t), (diff(y(t), t))(t), (diff(y(t), t, t))(t), (sin(y(t)+phi(theta(t))))(t), ((D@@2)(phi))(theta(t)), ((D(phi))(theta(t)))(t), (((D@@2)(phi))(theta(t)))(t).

Simple example: What does (diff(y(t), t))(t) mean?

g0;
IDTS:=indets(g0,function);
for ttt in IDTS do
  try
    res:=diff(ttt,theta(t))
  catch:
    print(lasterror);
  end try;
  Diff(ttt,theta(t))=res
end do;

for ttt in IDTS do
  try
    res:=diff(ttt,diff(theta(t),t))
  catch:
    print(lasterror);
  end try;
  Diff(ttt,diff(theta(t),t))=res
end do;

###################
This reminds me of an error some of my students would make occasionally, but which didn't seem to have consequences when plotting:

restart;
f:=sin(t); #f is not defined as a procedure (function) so should not be used as such.
plot(f,t=0..Pi); #Correct syntax
plot(f(t),t=0..Pi); # Incorrect, but works!
## So why does it work? Because numbers are accepted as constant functions:
5(7); #output 5
f(6);
subs(t=5,f(t));
%;

It seems that you want to find the jacobian matrix for the right hand sides of the system ode1,ode2 below.
I would do it like this:

restart;
f1:=(x,y)->x^2+y^2;
f2:=(x,y)->x^2-y^2;
ode1:=diff(x(t),t)=f1(x(t),y(t));
ode2:=diff(y(t),t)=f2(x(t),y(t));
VectorCalculus:-Jacobian([f1(x,y),f2(x,y)],[x,y]);

@roya2001 To determine 2 parameters you need 2 extra boundary conditions, not just one.

Let the loop round through this set K of 2 extra:

K:=combinat:-choose(extra_bcs, 2);

for b in K do
 try
   print(b=~0);
   res[b]:=dsolve(dsys4 union (b=~0), numeric, and the rest as is);
   ....
end do;

I got 11 successes.

My first reaction is that a thing like

(D@@2)(phi)(theta(t))(t); #(I removed superfluous parentheses)

is somewhat weird.

It means the second derivative of phi taken at theta(t) and then evaluated at t. That must be rethought!

1. You call the expression exp(-h^2/T)*(Int(exp(-x^2/T)*BesselI(0, h*x/T)*x, x = 0 .. 1))/T an equation. So is that expression supposed to be equal to zero?

2. What is the unknown: h or T?

4. BesselI(0,y) appears to be positive for all real y, so it seems you have a problem.

So what do you expect us to do?