fnormal(..., eps) and simplify(..,zero) as in the attached.

 restart:
  v1 := 1:
  v2 := 0.5:
  k := 0.12:
  p1 := 0.328:
  p2 := 0.74:
  l1 := 0.31:
  l2 := 1.16:
  sysode := diff(Cp(t), t) = -(p1 + p2 + k)*Cp(t) + l1*C1(t) + l2*C2(t),
            diff(C1(t), t) = (p1*Cp(t) - l1*C1(t))/v1,
            diff(C2(t), t) = (p2*Cp(t) - l2*C2(t))/v2:
  ics := Cp(0) = 0.347, C1(0) = 0., C2(0) = 0.:
  sol := evalf( dsolve([sysode, ics])):
  ans:=simplify(fnormal~(sol, 9), zero);
  plot( [ eval(Cp(t), ans), eval(C1(t), ans), eval(C2(t), ans)], t=0..10, color=[red, green, blue]);

Or, just do the whole thing numerically with

#
# Or do the whole thing numerically
#
  sol:=dsolve([sysode, ics], numeric);
  plots:-odeplot(sol, [ [t, Cp(t)], [t, C1(t)], [t, C2(t)]], t=0..10, color=[red, green, blue]); 

Download odesol.mw

to uses Newton's root-finding method?  ie https://en.wikipedia.org/wiki/Newton%27s_method

I have implemented this method for your problem in the attached

Download newton.mw

Rather than mahe a lot of assumptions, it is generally a better idea to use the evalc() command which evaluates complex expressions on th basis that all names represent real quantities, and does its very best to get a result in the form expr1+I*expr2

The square root of -1 is represented as 'I' in Maple , not 'j' (although you can change this if you really want to. In the attached I'm using the default.

It can be confusing to load a number of packages which are never actually used

Having said all of the above, the attached shows how I would perform this calculation.

Download refload.mw

you can only have four boundary/initial conditions - so in the attached I have commented out a couple of the ones in your original worksheet. No justification for the ones I commented out, except that one has to be on x(s) and one has to be on y(s).

After that (in order to obtain a numerical solution), numerical values have to be giocven to the names 'r' and 't' in your equations. I have (arbitrarliy) set these both equal to 1. Pick any other values you might want

With these changes a numerical solution i easily obtianed - see the attached

(And somebody has broken the inline display of worksheets on this site AGAIN)

Download odesol.mw
 

just by *looking* at the equations, that no solution is possible, which is what Maple returns. See the attached

Download noSol.mw

your are asking for a tangent plane at the point t=11, u=-12, it is probably a good idea to ensure that your plot ranges include this point.

The exponential function will be changing very rapidly at this point so you probably want to keeep the plotting ranges quite small - so something like the attached

Download tplane.mw

to passing the variable name as a parameter is to use the indets() function to extract the appropriate name from the passed polynomial as in the attached

Download quot.mw

functions with no errors - so what exactly is the problem?

You might want to check the execution group I have added at the end, but I doubt that this is what you are after

Download compRoots.mw

see the two attached worksheets below. Respectfully suggest that you submit an SCR

Download RG2020.mw

Download RG2022.mw

I mean what wiould you expect Maple to do with something which "looked" like the result you desire - see the attached for how to produce it - but what do you think you can do with it?

Download oddSeries.mw

'D' is not the differential operator - syntax is wrong. I assume Maple is just interpreting it as some (unknown?) function D(). The simple way to check is just to use convert(expression, diff) as in the attached, where I have also shown the correct syntax, and convert(expression, diff) works as expected.

This is all a bit "academic" because I don't think there is any waay to make the differetnial operator 'D' behave as if it were a "coefficient  of a polynomial'

Download Ddiff.mw

whether you are after a a "single" curve or multiple curves on the same graph. The attached does the latter

Download ptplt2.mw

whihc is perhaps closer to your original code, (and may be easier to understand?), is shown in the attached. All I really did was

1. Dispense with things that "looked" like functions, but weren't and really didn't need to be in the first place - so v(x), R(x), phi[n](x) etc can all be simply replaced by v, R and phi[n]. This avoids some notational confusion (for me at least)
2. Reordered the way the calculation is done by shifting the execution groups around into what I consider a more "logical" sequence (although your brain may work differently from mine!)
3. Never use sum() when add() will suffice.

Download oddSolve.mw

the function in the attached?

Download pw.mw

a case where the provision of a worksheet (use the big green up-arrow in the Mapleprimes toolbar) would be useful. Without such a worksheet, only general guidelines - and sort-of pseudoCode - can be given. Suggestions -

1. I'd split the problem into two stages, so that the "inner" maximization can be easliy tested/debugged
2. I'd force an optimization method which avoids the use of derivatives
3. I'd supply as much information in the form of constraints and/or variable ranges as possible

So the pseudoCode below would seem to be to logical approach to the problem. (NB this does not excute - it is pseudoCode)

  myProc:= proc( p1, p2, p3)
                 local Gamma, f, z1, z2, z3:
               #
               # Define the function Gamma
               #
                 Gamma:=someExpressionInvolving( f, z1, z2, z3);
               #
               # Evaluate this function for the 'passed' values
               #
                 myfunc:= eval(Gamma, [z1=p1, z2=p2, z3=p3]);
               #
               # myfunc is now a univariate function of the variable 'f'
               # It *may* be complex and may or may not be differentiable
               # so force an optimization method which doesn't use
               # derivatives
               #
                 sol:=NLPSolve( abs(myfunc),
                                f=f1..f2,
                                maximize=true,
                                method=nonlinearsimplex
                              );
               #
               # The above (if successful) will return a list containg
               # the required maximum, and the value of 'f' where this
               # was obtained. Ensure that this procedure returns only
               # the former
               #
                 return sol[1];
            end proc:
#
# Now one can test whether the above procedure is "functional", for a
# for a few sets of variables
#
  myProc(1.0, 2.0, 3.0);
  myproc(2.5e-6, 123.45, 2.0);
#
# Assuming the above test work (after debugging myProc if necessary),
# one can move on to the second stage of optimization, whihc would be
# something like
#
  NLPsolve( myProc,
            initialpoint=[1,1,1],
            method=nonlinearsimplex
          );
#
# If any kind of bounds on the variables are known, these should be
# supplied as well