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

@acer You are right (of course).
I just checked Fede's worksheet without making any changes and using Maple 18.02 on my old Windows 7, 32 bit machine. The timings were roughly the same. I encountered no problems at all with NumSteps= 6, 20, 40, 80.

It would help if you uploaded a worksheet. Use the green arrow in the MaplePrimes editor.

@acer Using method = _d01amc sure speeds up the computation: In Maple 18.02 (with the matrix<->Matrix replacement I mentioned) I got 0.44 seconds on my old machine and 0.17 on my new machine with Maple2015.2.

You use matrix instead of Matrix in the definition of A. That may make a difference. (Apparently not: see below).
I changed to:

A := Matrix(NumSteps, 2, datatype = float);

Similarly, I added the datatype option to the definitions of T1 and AF.

I measured the timing of the loop only.
For NumSteps=40 in Maple 18.02 on an old not very fast machine I got 21 seconds. For NumSteps=80 the time used was 44 seconds.

In Maple 2015.2 on a new and reasonably good machine the times were about half these amounts.

@want to be a permanent vegan You refer to 3 worksheets. I had a brief look at one of them (example_hw2.mw) , which has these integrals:

and

Since p(t) apparently doesn't depend on x these two integrals are obviously easily computed, just as in

restart;
int(p(t),x=-eta(t)..eta(t));
                             

Surely you don't intend anything as trivial as that?
So I'll stick to my statement about mathematics: Formulate the problem mathematically, so it is obvious to any person with some background in mathematics.
The point is (as somebody said): If you don't know where you want to go, any road will take you there.

@crazyeti You wrote

"I just hope Maple could post a solution/workaround soon, something like always executing the assignment statement twice etc."

That won't work as it is the second one that is wrong.

Added: I have no idea about what is going on, but it is somewhat entertaining (Maple2015.2):

restart;
expr:='a+b+c+d+e+f+g+h+i+j+k+l+m+y-y+n+o+p+q+r+s+t+u+v+w+x+z';
expr;
expr;
indets(expr,name);
nops(%);
nops(expr);
eval(expr,{w=W,z=Z}); #2*W and 2*Z