You could call rsolve on the system, after converting floats to exact rationals. This produces an explicit solution.

You could then utilize the explicit result, eg. applying evalf to the result formulas, or applying evalf after using values for t.  See each of ansE, andF, and ansFF below. Of those, offhand I would suggest the first scheme with ansE, substituting values of t into the exact formulas and then applying evalf as a last step. Your mileage may vary.

You could experiment with higher Digits, when evaluating and evalf'ing at concrete values of t.

Download SystemRecursive_ac.mw

You could compare results with those from a recursive procedure such as Kitonum supplied. Don't expect both to agree to all decimal places.

The following attachment demonstrates re-using the "preconditioned" information (eg. column switching to minimize fill-in, etc), given that the new Matrix has the exact same sparsity pattern as the original.

sparse_direct_symb8.mw

The next attachment replaces my SLU procedure from my previous code, as an edit to your revision of the original worksheet. It implements the previous SLU's actions but in a finer set of steps, allowing the "symbolic precontioninging" to be reused for subsequent calls.

Phasefield2DBaseCodeParametricMay02_ac1.mw

It gets some timing improvement over my previous. But I have not had time to replace the Jac1 procedure with an efficient scheme of generating the sparse matrix's triple-of-vectors (neither compressed nor uncompressed row-column index vectors). As the problem size grew that Jac1 incurred a significant portion of the time cost.

As bonus, UMFPACK's Control array is here directly available, and you could even try tweaking certain parameters of the solver. (Earlier I had mistakenly presumed that array was a struct, so it was easier than I'd imagined.)

Perhaps you could utilize this (slightly modified, accordingly) alongside your code that efficiently produces the sparse Matrix directly in a triple-vector format.

[edit] For anyone else reading this, the goal here is to remove the need for calling umfpack_dl_symbolic following the first sparse linear solve. This is remove duplicated effort in constructing the preconditioning -- given that the sparsity pattern of the lhs matrix remains the same for all linear solves.

The goal here is not to provide functionality to allow different rhs B in repeated linear solving. In fact since I didn't provide a way to undefine the module's Sym and Num locals then at present it can only be used for a single sparsity pattern per restart. And it also doesn't have functionlity to "free" memory associate with those handles.

You can use the eval command for this, to evaluate expressions using substitutions for x.

These eval calls can each be passed directly as the first argument in calls to the plot command.

V1 := -x:
V2:= -2*x + 4.95:
V3:= -x + 1.665:

l_1:=1.665: l_2:=4.95:

eval(V2, x=x-l_1);

        -2 x + 8.280

eval(V3, x=x-l_1-l_2);

         -x + 8.280

There are several variants on all that. Here is one such, in which l_1 and l_2 do not need to actually be assigned.

restart;
V1 := -x:
V2:= -2*x + 4.95:
V3:= -x + 1.665:

params := [ l_1=1.665, l_2=4.95 ]:

eval(V2, x=eval(x-l_1, params));

        -2 x + 8.280

eval(V3, x=eval(x-l_1-l_2, params));

         -x + 8.280

Download add_map_sets.mw

The following uses RootFinding:-NextZero (internally):

   fsolve(cos(x)=x^2/1000,x=-40..40,maxsols=infinity);

On my Maple 2022.0 (and 2018.2, which is what the OP seems to have used) that returns 22 roots.

This is a common mistake.

One way to deal with it is to construct the operator (procedure) G by using the unapply command.

Another way is to use expressions instead of operators for g and G.

Download expr_vs_operators.mw

f:=(x,y)->piecewise(x<>0 or y<>0,(y^2+x^2)^x,1):

DX,DY:=D[1](f),D[2](f):

DX(1,1),DY(1,1);

               2 ln(2) + 2, 2

If you are looking to construct the following (or similar) then you could utilize the foldl command instead of rolling your own implementation.

L:=[a,b,c,d]:
v:=[p,q,r,s]:

foldl((z,j)->F(z,L[j])-v[j],y,$1..nops(L));

   F(F(F(F(y, a) - p, b) - q, c) - r, d) - s

I added option redraw=false to the plots:-display calls, and option adaptive=true to the plot calls.

I also changed sum to add, and Int to Int, and a few other minor edits.

The result runs in my 2022.0 about 5% faster than in 2020.2.

marlin01.mw

Here's one way:

Download rsolve_fun.mw

I did that in Maple 2022.0. In, say, Maple 2020.2 I could wrap evalc(simplify(ansX)) to get the pure trig form instead of terms like I^k and (-I)^k. I guess that aspect of the form is a matter of preference.

Download sum_simp.mw

That temp1 could also be attained by the following.

subsindets(simplify(ee),specfunc({sum,Sum}),
           u->expand(op(0,u)(simplify(content(op(1,u)))
                             *simplify(op(1,u)/content(op(1,u))),
                             op(2..,u))));

The later steps are just for formatting of an inert multiplication by -5/3.

Other examples might require slightly different formulations -- I can't tell what might be your more general case (... if there is such).

Here are two similar approaches, both using `%*` for an inert multiplication. The second gets pretty-printed nicely under both extended and standard typesetting level, whereas the first does so only for extended.

Download yank.mw

That numeric (rational) factor will get distributed across that sum (when that N is not also put into the denominator). That happens due to automatic simplification, which occurs between parsing and evaluation and which thus cannot be prevented by unevaluation quotes.

'4*(x-y)';

                4 x - 4 y

In the solution above I've prevented the automatic distribution by using an inert multiplication. There are other ways around it. One is to wrap the numeric factor in a call to ``(...), ie. ``(1/1296) but I don't like thay because it renders wrapped in extraneous round brackets.

lprint(indets(sys, name));

   {x, c[0], c[1], c[2], c[3], eqq[0], eqq[3]}

Your last loop only assigns to eqq[1] and eqq[2]. Your code has that loop as,

for i to n-1 do
  eq[i] := RHS_weak_form(i);
end do

Perhaps you intended something more like the following (I am guessing):

for i from 0 to n do
  eqq[i] := evalf(eq[i]-RHS(i));
end do;

But if that is the case then your system would have five unknowns: {x, c[0], c[1], c[2], c[3]}, but x is a dummy name of integration. By informing fsolve that the variables are only {c[0], c[1], c[2], c[3]}, I get the following,

SOL := fsolve(sys, {c[0], c[1], c[2], c[3]});

    {c[0] = 1.000000000, c[1] = -13.06296523,
     c[2] = -11.79236213, c[3] = 2.718281828}

Revised worksheet:

Download problem_solve_system_ac.mw

Either rationalize or evalc can get that task done. And then you might further adjust as you wish.

Download simpex.mw

note: You used rationalize to similar effect in a previous Question.

Here is another way, with a single call to solve.

Its result fortuitously contains the pair of solutions with calls to 2-argument arctan, as you wanted.

Each also happens to contain the beta(t) solution, but if needed then it's trivial to sieve that out -- as shown in the last line below.

Download Arctan_ac.mw

You're original example happened to get the intermediate solution (for beta(t)) unsimplified, because subs does not induce a final evaluation.

Your undesired example was not due to automatic simplification, which has a very specific meaning. It was due to the usual full evaluation that happens to arguments of procedure calls, eg. calls to solve. You can also pass the intermediate result from subs -- referenced as its own separate result -- in an argument of a call to solve, while avoiding the usual full evaluation by means forced 1-level eval.  Arctan_ex.mw

But this is all quite problem-specific. I don't see how you can reasonably expect always to be able to get semi-generic solutions that just happen to cover (explicitly, separately) only a pair of branches.