You were failing to provide any knowledge of f to the attempts at computing c5,c6,c7.

I resolved that issue by calling dsolve on a coupled system -- once, for each set of parameter values. (Yes, you could dsolve for f alone, decoupled, and then reuse that result in the later dsolve calls for theta. But efficiency seems not a major concern in what you've attempted here. I'm just trying to keep it simpler for you to do -- it's just a judgement call, on what might be best for you.)

In at least one place you also had f instead of f(x). That's the kind of mistake you can make when you allow declare to make things visually unclear.

You had duplicates of most of the parameter information. I've cleaned all that up and made it more organized, which might help you keep track of which is what. Having twice as many solutions as necessary floating around, with numbering mismatches, seems a recipe for confusion.

Note that M only appears in eq1, and is multiplied by sin(alpha).
Hence M has no effect when sin(alpha) is taken as 0, which is the case for all your given parameter sets.

So I rearranged the parameter sets according to the value for cos(omega). And I added a legend. We can see that for your parameter sets it is the value of cos(omega) that is the key difference.

Download IP-TEMP_ac.mw

You might wish to also try varying the sin(alpha) parameter, and then see the related effects of varying M.

The keen observer might notice that there is a lot in common amongst all those dsolve and odeplot calls. There is so much that could be generalized. But I've deliberately chosen not to wrap it all up in a single resuable procedure that could produce all the results.

There are a lot of members on this forum who attempt very similar problems. Several struggle with numeric dsolve & odeplot, and perhaps even with the DE concepts. I think that some members can learn more by doing each step themselves rather than use an abstracted do-it-all solution procedure.

The following was done using Maple 2023.0.

Where you see R subscripted by x it was entered as the keystrokes R__x with two underscores, and similarly for y and z subscripts.

Download sincos_ts.mw

The code in red above could be tucked away out of sight in a worksheet's Startup Code, accessible via the main menubar's Edit menu.

ps. I could mention that the first three lines of red code, which assigns to (s,c) and the print-extensions, could be replaced by alias, as Edgardo mentioned. My longstanding preference is not to use alias, which is why I didn't mention it earlier. The rest of the red code, the Typesetting calls, accomplish the rest of the aspects you mentioned or showed.

That lengthy result you obtained appears to have a "hidden" division-by-zero in it. I suspect that you might not be able to just simplify that away, given that form.

But the option method=trager to the int command seems to help here.

Download int_tr_ex.mw

I don't know all the details of the control flow by which available methods are attempted by default for indefinite integration. I'll submit this example as a bug report.

Consider this counterexample, where x=-1.

Download cbrt.mw

Perhaps a suitable time to review the documentation for the root (esp. 4th and 5th bullet points) and surd commands.

This particular example can be handled in one call to fsolve.

Download rf_ex.mw

note: Given a numeric rootfinder, some example can be found that defeats it.

The contents of the Favorites palette are stored in a location on one's Operating System, somewhat similarly to how one's Maple GUI preferences are saved.

They are not saved inside an individual .mw file.

Also, Maple preferences are stored on the local computer, say, not in the cloud.

The contents of the palettes is part of the GUI itself, not the document.

Once you've created the plot using geometry:-draw, you could merge it with other kinds of plots in the usual way, ie. using plots:-display.

Download plotpoint_ac.mw

You can replace those function calls (of t). You could do it using evalindets, and also with eval and a set/list of replacements.

This targets those particular names, but it'd also work if you had terms like sin(t), etc, leaving those alone.

problem_ach_ac.mw

Download vec_dot_add.mw

Is this the kind of thing you meant?

Download comp_map.mw

Do you mean something like this, to get that data into the columns of B?

  M0:=Edges~(graphs);
  M1:=Girth~(graphs);
  M2:=Diameter~(graphs);
  B := <<M0>|<M1>|<M2>>:

And then, as you had it,

  Export(..., B);

or, say,

  Export(..., B, format="delimited"):

or some flavour of call to ExportMatrix, etc.

Here are three straightforward ways to map that operation across your ptlist,

   (u->TransMatrix.u)~(ptlist);

   map(u->TransMatrix.u, ptlist);

   map[2](`.`, TransMatrix, ptlist);

Are you trying to accomplish something like this, using the command Physics:-diff to differentiate wrt the function call u(t)?

note: what you showed in plaintext is not valid 1D Maple Notation code, so I've guessed at what you might have intended. Adjust as needed.

restart;
H:=(y,u,t)->y(t)^2+u(t)^2+1/(2)(x(t))^(4) +(`∂/∂ x`)(V(t))*u(t):

Physics:-diff(H(y,u,t), u(t));

   2*u(t)+(`∂/∂ x`)(V(t))

H:=(y,u,t)->y(t)^2+u(t)^2+1/(2)(x(t))^(4) + Diff(V(t),x)*u(t):

Physics:-diff(H(y,u,t), u(t));

        2*u(t)+Diff(V(t),x)

You haven't shown us your actual expression. But you might be able to temporarily substitute for x-x0, then collect, then resubstitute.

Btw, your attempt using abs(x-x0) works as a specialized case of collecting wrt just abs.

Download collect_algsubs.mw

The results from implicitplot3d often look clunky, and sometimes edges inaccurately appear rounded, unless the sampling grid option is set so high that the GUI becomes sluggish when they are manually rotated, etc. It is an inefficient mechanism. It is a poor choice for handling this example.

I prefer to use plot3d when possible. Or, in this case, the tetrahedron plotting command can be used and produce an even more lightweight result.

Both of those can produce a much better rendition here than does implicitplot3d. 

Download tetra.mw

Unfortunately shadebetween cannot handle a variable range end-point. Hence my alternative was to fill a plot3d of y-1 and then raise by 1.