@acer In the meantime, if there is a strong enough need for such a facility, then one could use IntegrationTools:-CollapseNested for this purpose:

CollapseNestDiff:=proc(u::specfunc({diff,Diff,%diff})) local a,b;
   a:=subs(diff=Int,Diff=Int,%diff=Int,u);
   b:=IntegrationTools:-CollapseNested(a);
   subs(Int=Diff,b)
end proc;

@John Fredsted That switching surprised me, but take a look at the lprint versions:

expr := Diff(f(x,y),x,y,y):
lprint(expr);
Diff(f(x, y), x, y, y)
expr1 := convert(expr,diff):
lprint(expr1);
diff(diff(diff(f(x, y), x), y), y)
expr2 := convert(expr1,Diff):
lprint(expr2);
Diff(Diff(Diff(f(x, y), x), y), y)

@Honigmelone You may also want to look at Physics:-diff:

restart;
a:=b(t)+c(t)*diff(b(t),t)+diff(b(t),t$2)+diff(b(t),t$3)+45*diff(c(t),t);
Physics:-diff(a,diff(b(t),t));
Physics:-diff(a,diff(c(t),t));
Physics:-diff(a,diff(b(t),t,t));

@iman What kind of symmetry are you thinking about? A graph of a function of x couldn't possibly be symmetric about the x-axis unless the function is zero on the whole interval.

The solution you are getting I wouldn't trust though.
Try doing:
SYS:=expand(solve(newsys,{diff(g1(y),y$3),diff(g2(y),y$4),diff(g3(y),y$4)}));
and observe the denominators on the right hand sides.

You can expect problems for y=0.
You set abserr=1e10, i.e. 10^10. That is kind of large for an absolute error where the values of g1,g2,g3 are in fact pretty small (~10^(-9)).

@iman L is never given a value. What is it?

@iman No in my two uses of dsolve above a finite difference method is used.
Shooting as I understand it works like this on a boundary value problem.
Again an example:
restart;
sys:={diff(y(x),x,x)+lambda(x)*y(x)=0,diff(lambda(x),x)};
bcs:={y(0)=0,y(Pi)=0,D(y)(0)=1};
##The finite difference result first:
res1:=dsolve(sys union bcs,numeric);
# Now a shooting method:
ics:={y(0)=0,D(y)(0)=1,lambda(0)=a}; #A parameter a to be determined
##Notice that ics are initial conditions, not boundary conditions.
## We must determine the parameter a such that y(Pi)=0:
res2:=dsolve(sys union ics,numeric, parameters=[a],output=listprocedure);
Y:=subs(res2,y(x)); #Procedure to evaluate y(x)
Q:=proc(a) res2(parameters=[a]); Y(Pi) end proc;
L:=fsolve(Q,0.1234); #even without the guess it works: fsolve(Q)
res2(parameters=[L]);
res2(0);
res1(0);

@iman
I believe that the method used by dsolve/numeric/bvp to find a parameter is the way I do it in the following simple example (the second way below). I know of no name for that method, but it is pretty straightforward.

restart;
ode:=diff(y(x),x,x)+lambda*y(x)=0;
res1:=dsolve({ode,y(0)=0,y(Pi)=0,D(y)(0)=1},numeric);
res1(0);
## Now the explicit way of doing the same:
sys:={subs(lambda=lambda(x),ode),diff(lambda(x),x)=0};
bcs:={y(0)=0,y(Pi)=0,D(y)(0)=1};
res2:=dsolve(sys union bcs,numeric);
res2(0);

@iman The method used used by dsolve/numeric/bvp is a finite difference method, not a shooting method.

Besides the problems mentioned by acer, it seems that no value is ever given to a, which appears in extra_bcs.

Switching to 1D input (as suggested by acer) will make it easier for me (and I believe you) to understand what is going on.

What was the intended meaning of this last  one of your set ODE:

diff(y(eta), eta)+Sb*f(eta)*y(eta)-Pe*(w(eta), eta)*(diff(t(eta), eta))-Pe*t(eta)*y(eta) = 0;

Instead of (w(eta), eta) did you want just w(eta) or did you want diff(w(eta), eta) ?

If I assume the former (i.e. w(eta) ) then by doing

etainf:=15:
BCS2:=BC union remove(hastype,IC,name);
dsolve(ODE union BCS2,numeric);
# I get the error message:
Error, (in dsolve/numeric/bvp/convertsys) too many boundary conditions: expected 12, got 15

@tomleslie I think the comment may have been meant for me, but put in a wrong place.
I tried to answer his question as he put it: "Problem with seq", i.e. as a problem having to do with seq.
Then he had the (to me) puzzling comment to my answer saying: "However, the output gives commas and quotes".
I now see from this latest comment that he is (now at least) using print and not printf.

@das1404 I have added an example below. Your comment "However, the output gives commas and quotes" puzzles me.

@tomleslie I never used that icon, but just tried entering sqrt(1-u) after no success with sqrt(1+u).

No success with sqrt(1+sqrt(5)), but change the sign and voila:

No success with sqrt(3+sqrt(5)), but sqrt(3-sqrt(5)) gives

I shall not bother you with more examples although it is amusing to try and see what works.

What is it that is being repeated? I see absolutely nothing being repeated.

@iman Actually, my comments were intended to dissuade you from pursuing the symbolic approach and take up the numerical approach instead.

@iman Before you get started on solving the problem numerically, be sure to replace omega^2 by omega2:
sys:=subs(omega^2=omega2,{p1,p2,p3});