which you can tell more or less by inspection of the form of the PDEs and the boundary conditiions

You can persuade Maple to provide these answers using the code below

  restart;
  PDE1 := diff(u(y, t), t)+diff(u(y, t), t, t) = diff(u(y, t), y, y)-u(y, t);
  PDE2 := v(y, t)+diff(v(y, t), t) = diff(u(y, t), y);
  ICandBC := {u(0, t) = 0, u(3, t) = 0, u(y, 0) = 0, v(y, 0) = 0, (D[2](u))(y, 0) = 0};
#
# Note that PDE1 is independent of v(y,t), so can be solved
# independently of PDE2 - and can in fact be solve analytically
# Unfortunateely it returns u(y,t)=0, which I suppose one should
# have guessed from the form of the PDE and the ICs/BCs
#
  solpde1:=pdsolve( [PDE1, ICandBC[1..3][], ICandBC[5]]);
#
# Substitute this solution in the second PDE
#
  PDE2:=expand(subs(solpde1, PDE2));
#
# Now solve this "new" second PDE, which Maple can't
# do analytically but can do numerically
#
# Unfortunately this too returns u(y,t)=0
#
  solpde2:=pdsolve( PDE2, [ICandBC[4]], numeric, range=0..3);
  solpde2:-plot3d( v(y,t), y=0..3, t=0..5);

The only part which seems relevant to "building a hyperbola" is the pharse " show that the place of C when a varies is a branch of hyberbola"

The following code draws the position ofyuor point  "C" as "a" varies

  restart;
  with(plots):
  alpha:= arccos((2*L-a)^2/(2*a*(2*L-a))):
  h:= tan(alpha)*(2*L-a)/2:
  C:= eval~([L-(2*L-a)/2, h], L=6);
#
# The following "looks" like a hyperbola
#
  display( [ plot( [ C[1], C[2], a=4..10]),
             plot( [ C[1], -C[2],a=4..10])
           ]
         );
#
# Or swapping x- and y-coordinates
#
  display( [ plot( [ C[2], C[1], a=4..10]),
             plot( [ -C[2], C[1], a=4..10])
           ]
         );

which produces

and the attached now produces plots - whether they are the plots you want?????

plt.mw

which produce the plots

Rarely have I seem such confusion betwee what is a 'function' and what is an 'assignment' or what is a variable and what is a 'parameter'. Makes figuring out what you actually want almost impossible. So the attached is just a wild guess - but at least it is organsised 'logically' so you ought to be able to refer to it when explaining what you actually do want!!
 

  restart:
#
# Define a load of parameters
#
  beta:=2:
  alpha:=2:
  sigma:=3:
  A:=3.3:
  p:=1:
  q:=1:
  rho:=4.5:
  r:=1:
  v:=4:
  Theta:=(beta^2)-4*alpha*sigma:
  xi:=x-(2*r/(log(A)^2*Theta*(v^2-1)))*y-v*t:
  tanA:=I*((p*(A^(I*xi))-q*(A^(-I*xi)))/((p*(A^(I*xi))+q*(A^(-I*xi))))):
  Q:=-beta/(2*sigma) + sqrt(4*alpha*sigma - beta^2)/(2*sigma)*tanA*(sqrt(4*alpha*sigma - beta^2)/2*xi):
#
# Define some equations
#
  EQ1:= c*beta*alpha*b[1]*(v - 1)*(v + 1)*ln(A)^2 - 2*v*b[0]^3 - r*b[0]:
  EQ2:= b[1]*(c*(v - 1)*(v + 1)*(2*alpha*sigma + beta^2)*ln(A)^2 - 6*v*b[0]^2 - r):
  EQ3:= 3*b[1]*(c*beta*sigma*(v - 1)*(v + 1)*ln(A)^2 - 2*v*b[0]*b[1]):
  EQ4:= 2*(c*sigma^2*(v - 1)*(v + 1)*ln(A)^2 - v*b[1]^2)*b[1]:
#
# Solve the system of equations. Note that there are five
# solutions. Only the fourth and fifth solutions give
# *explicit* values for all of the unknowns
#
  sols:=solve({EQ1, EQ2, EQ3, EQ4},{b[0],b[1],c},explicit);
#
# Plot the fourth and fifth solutions of the above - ie
# the only two which do not involve arbitrary values for
# some parameters
#
# Note that it is necessary to ovide a value for the name
# 't'. Should this quantity have been defined in the 
# parameter list above? Who knows?
#
  plot3d( eval(eval( b[0] + b[1]*Q, sols[4] ), t=1),
           x=-10..10, y=0..1, size=[1000, 1000]
        );
  plot3d( eval(eval( b[0] + b[1]*Q, sols[5] ), t=1),
          x=-10..10, y=0..1, size=[1000, 1000]
        ); 

which produces the figures

and

See the attached

sys.mw

If you don't believe me, try

plot(P1, x = -10 .. 10)

(neither better nor worse), but just because I like alternatives

restart;
  M:=Matrix(3, 3, [ [-1/3*a+b+c-2*d, -1/4*a+2/3*c, 1/15*a-1/4*c+2/3*d ],
                    [    e+a-2*c,        2/3*a,       -1/4*a+2/3*c    ],
                    [     f-2*a,           0,             2/3*a       ]
                  ]
           ):
  Id:=LinearAlgebra:-IdentityMatrix(3):
  solve([entries(M=~Id, 'nolist')]);

but if I fix some issues related to variable scope, the attached

pts.mw

will produce the following plot

When you capitalise the sum() function, you are making it "inert". The simple change in the code below provides an answer

 restart;
 rho:=1/(2);
 mu:=1/(2);
 T:=1/(4)*(t-x)^(1/(2));
 E[rho,mu](T):=sum((T^(k))/(GAMMA(k*rho+ mu)),k=0..5);
 eq1 := h(t) = (8/3)*(int((t-x)^(1/2)*E[rho, mu](T)*(1/(sqrt(Pi)*x^(1/2))-(3/32)*x^2+(3/16)*exp(-x)-1), t = -1 .. 1));

Working correctly the attached

dsol.mw

because you posted a 'picture' of it which is unusable (and I'm not going to retype it.

However you may get something from considering the toy example in the code and attachment below

  restart;
  with(Statistics):
  A := < seq(sin(2*Pi*i/50), i = 1..15) >:
  B := < seq(0.5, i = 1..15) >:
  C := < seq(0.05, i = 1..15) >:
  ErrorPlot(A, xerrors = B, yerrors = C);

which will produce the graph below (although there are many options, bells and whistles which you can include in the ErrorPlot() command

scat.mw

is show below (Actually pretty similar to Kitonum's solution!)

restart;
  L0:=[w,b,g]:
  L:=combinat:-choose([seq(j$4, j in L0)], 4);
  occ:=(el->[seq(s=ListTools:-Occurrences(s,el), s=L0)])~(L);

which produces

 L := [[b, b, b, b], [b, b, b, g], [b, b, b, w], [b, b, g, g], 

  [b, b, g, w], [b, b, w, w], [b, g, g, g], [b, g, g, w], 

   [b, g, w, w], [b, w, w, w], [g, g, g, g], [g, g, g, w], 

   [g, g, w, w], [g, w, w, w], [w, w, w, w]]


     occ := [[w = 0, b = 4, g = 0], [w = 0, b = 3, g = 1], 

       [w = 1, b = 3, g = 0], [w = 0, b = 2, g = 2], 

       [w = 1, b = 2, g = 1], [w = 2, b = 2, g = 0], 

       [w = 0, b = 1, g = 3], [w = 1, b = 1, g = 2], 

       [w = 2, b = 1, g = 1], [w = 3, b = 1, g = 0], 

       [w = 0, b = 0, g = 4], [w = 1, b = 0, g = 3], 

       [w = 2, b = 0, g = 2], [w = 3, b = 0, g = 1], 

       [w = 4, b = 0, g = 0]]

Maple will solve your boundary-value problem directly, so converting it to an intial-value problem (the essence of the shooting method) seems rather pointless!

Such conversion from BVPs to IVPs is gnerally only used for ODE systems whic are being solved numerically. (As a general rule numerical IVP methods are "more robust" than numerical BVP methods)

The attached solves your original boundary-value problem, then converts it to to an initial-value problem and solves the latter. Same answer both ways

  restart;
#
# Solve the original system as a boundary-value problem
# ie using the boundary conditions
#
# u(sigma) = 0, u(-sigma) = h
#
  ode := diff((diff(u(y), y) + epsilon*diff(u(y), y)^3)/(1 + epsilon*diff(u(y), y)^2), y) = 0;
  bcs := u(sigma) = 0, u(-sigma) = h;
  sol0 := dsolve([ode, bcs], u(y));
#
# Solve the same system using only one boundary condition,
# namely,
#
# u(-sigma) = h
#
# Note that the solution returned will contain an arbitrary
# constant
#
  bc1:= u(-sigma) = h:
  sol1:= dsolve([ode, bc1], u(y)):
#
# 1. Compute the value of the arbitrary constant which ensures
#    that the second boundary condition above is satisfied
# 2. Differentiate the solution, and evaluate it at the given
#    boundary condition ie u(sigma) = 0, using the value of
#    of the arbitrary constant obtained in (1) above.
#
#  This will generate a further boundary condition of the form
#
#  D(u)(-sigma) = someValue
#
  bc2:= convert(subs(isolate(rhs(eval(sol1, y = sigma)) = 0, _C1), eval(diff(sol1, y), y = -sigma)), D):
#
# Solve the same system with the boundary conditions
#
#  u(-sigma) = 0
#  D(u)(-sigma) = someValue
#
# Note that since both boundary conditions are stated at the
# same value of the independent variable, the ODE system
# is now na "initial value" problem, rather than a "boundary
# value" problem.
#
# Such conversion of BVPs to IVPs, is the essence of the "shooting"
# method
#
# The "new" boundary condtions for the IVP are
#
  bc1;
  bc2;
#
# Solving the same ode with these new boundary conditions
# generates the same answer as before
#
  sol2:= dsolve([ode, bc1, bc2], u(y));
#
# Check
#
  is(sol0=sol2);

shoot.mw

which really should not be necessary!!

Your problem does look like a *bug* in the UnderlyingGraph() command

  restart;
  with(GraphTheory):
  G1:= Digraph( [a,b,c],
                {[[a,b],2],[[b,c],3],[[c,a],4]}
              );
  Under:= proc(g::Graph)
               uses GraphTheory:
               return Graph( Vertices(g),
                             { seq
                               ( [ convert(j[1], set), j[2] ],
                                 j in Edges(g, weights=true)[]
                               )
                             }
                           ):
          end proc:
  G2:=Under(G1);
  Edges(G1, weights=true);
  Edges(G2, weights=true);

the code below

restart;
sol:=solve(0 < M^2 - 4*m^2, M) assuming (0 < M, 0 < m);
plots:-inequal([sol[1][],sol[2][]], M=-10..10, m=-10..10)

whihc produces the plot

slightly differetn way is shown in the attached

ellip.mw