@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});

In addition to the suggestion from Markiyan notice that odeplot allows functions of the dependent variables too, like ln(y(t)).

res1:=dsolve({op(eval(DissMod1,Pars1))} union {y(0) = 1, B[1](0) = 1, B[2](0) = 0}, numeric);
plots:-odeplot(res1,[t,y(t)], 0..1000);
plots:-odeplot(res1,[t,ln(y(t))], 0..1000);

@Carl Love That is indeed a good question.
I tried these:
restart;
ode1:=diff(x(t),t)=1/[x(t)+1]*sin(t);
res:=dsolve({ode1,x(0)=0},numeric); #No apparent problem
plots:-odeplot(res,[t,x(t)],0..1); #Now the problem shows up
## Also of course it would show up just by doing
res(0.1);
Error, (in res) cannot evaluate the solution past the initial point, problem may be complex, initially singular or improperly set up

ode2:=diff(x(t),t)=[x(t)+1]*sin(t);
dsolve({ode2,x(0)=0},numeric); #Problem right away
## It is not surprising that these two complain right away:
dsolve({ode1,x(0)=0},numeric,range=0..1,output=piecewise);
dsolve({ode1,x(0)=0},numeric,output=Array([0.5,1]));
##
My guess would be that dsolve/numeric (with output=procedurelist/listprocedure/operator) sets up a procedure (or several) to calculate the solution, but does the minimal amount of calculation to do so. Thus it is left to chance whether an expression like [x(t)+1] (which is not syntactically wrong) is caught.

Addendum.
Setting infolevel[dsolve]:=5 shows that the computation procedure for the ode1 problem has YP[1] = sin(t)/[Y[1]+1]  (coming from DEtools[convertsys]).
That problems come up right away for ode2 is due to the result of eval(rhs(ode2),t=0), which is [0].
For ode1 we get eval(rhs(ode1),t=0), which is just plain 0.
A third version ode3:=diff(x(t),t)={x(t)+1}*sin(t);
This one also survives dsolve because eval(rhs(ode3),t=0) is also just 0.