@Earl The plot in your worksheet is much too big for my computer to handle. The GUI alone consumed 13 Gig of memory, and I needed to kill my entire Maple session.

However, I can see from your code that you are using a coloring technique completely different from what I described.

@acer I agree in general with all that you said. However, you might not realize that my only purpose (and I'm guessing VV's only purpose) was to read (as a human) the code of partition1 out of mathematical interest. I wasn't trying to extract the code to execute it. I think that both of our methods are useful for just reading the code.

Whenever I'm on a code-reading exploration like that, I habitually set verboseproc to 2 or 3 and opaquemodules to false (whether it ends up being needed or not).

My comment about browsing the library was meant to be a remark on the difficulties one faces when trying to read the library (just as a human reader, not for execution purposes). I didn't mean to suggest that some large percentage of those numeric-name-only entries were similar to partition1.

By the way, why is the library mentioned specifically in that bit of the Inert tree that I extracted?

I guess that you're avoiding piecewise because you can use neither it nor Heaviside in off-the-shelf Fortran?

@vv Here's another way to get to the code of that procedure:

interface(verboseproc= 3):
kernelopts(opaquemodules= false): #not strictly needed in this case
indets(
    ToInert(eval(combinat:-partition)),
    And(
        specfunc(_Inert_ASSIGNEDLOCALNAME), 
        satisfies(x-> op(1,x)="partition1")
    )
);
{_Inert_ASSIGNEDLOCALNAME("partition1", "PROC", 
  36893490238314748452, _Inert_ATTRIBUTE(_Inert_STRING("C:\Progr\
  am Files\Maple 2021\lib\maple.mla")))}

P:= pointto(op([1,3], %));
                        P := partition1

eval(P);

Browsing through the main library with LibraryTools, I chagrinfully note that the vast majority (about 90%) of entries are anonymous, with their only "name" being an integer of 5 or 6 digits.

@vv If you replace option autocompile with option remember, then for the largest of the cases that you showed, the *non*-compiled procedure is 169 times faster than the compiled procedure.

@tomleslie I've seen that problem also in Maple 2021. For me, it didn't occur with this worksheet. I suspect that this bug has nothing to do with one's code and instead is based on environment 

@Preben Alsholm As far as I can see, it works in Maple 2021. Which part isn't working for you?

@vsnyder You wrote:

• Another problem (that isn't mentioned in the Maple Advanced Programming Guide) is that sending statements to EvalMapleStatement isn't quite the same as typing them into the Maple textual interrace. In the Maple textual interface, you can lay out a "proc" definition on several lines, as you would in C or Fortran or Haskell. To submit it to EvalMapleStatement, it needs to be combined into one long stream of text, from "proc" to "end proc".

That is completely false. You either haven't understood anything that I've said in this thread regarding line breaks, white space, etc., or you've chosen to ignore it.

@MapleEnthusiast You have a system of 4 equations in 4 unknowns, so I don't know why you think that its solution has a free variable that you can arbitrarily set to 1. 

It is much, much, much easier to find numeric solutions of sysP with fsolve than it is to work with the huge symbolic solution I generated in the Answer. For example:

fsolve(sysP, {vz[]}=~ 0.1..infinity);
        {z[2, 1] = 1.204997696, z[2, 2] = 1.168035195,  
         z[4, 1] = 3.038404729, z[4, 2] = 3.170903246}

Remember that these need to be squared to get the corresponding y values.

@elonzyy You wrote:

• to_list := a -> convert(a, list);
  

Better (for generic a): 

to_list:= `convert/list`:

If you're expecting a to be a certain type (Vector, table, set, etc.), then there are better ways.

• map_i := (f, a) -> zip(f, to_list(a), [`$`(1..nops(to_list(a)))]);

Yes, that's bloated due to unnecessary quotes and parentheses and inefficient due to the repetition of to_list. Better:

map_i:= (f,a)-> (L-> f~(L, [$1..nops(L)]))(to_list(a)):

• flat_map_i := (f, a) -> ListTools:-Flatten(map_i(f, a), 1);
  flat_map := (f, a) -> ListTools:-Flatten(map(f, a), 1);

Those two aren't even worth typing or giving a name to. You use an f that creates a sublist from a sequence, and then you turn that sublist back into a sequence. It'd be better to simply not create any list with f in the first place, for example,

map_i(1@@0, [3,5,7]);

(1@@0 is the multi-argument identity function. An equivalent form for it is ()-> args.)

@Zeineb The reason that your instructor assigned you this exercise is specifically to teach you some of the same things that I've been telling you. This function was specifically designed to cause trouble (an erratic sequence) for Newton's method when the initial value is several "humps" away from the true solution. Look at a plot of f and another of its derivative on the interval x= 0..22. Note the wild oscillations of the derivative. If the derivative is near zero, the tangent line is nearly horizontal, and thus the point where the tangent intersects the x-axis (which is, by definition, the next x) will be far away from the current x.

Here's a procedure that I think will help you see the "bare bones" of what's happening. You can use it as an alternative to Newton in Student:-NumericalAnalysis. You should especially experiment with different values of digits and different initial values.

Newton:= proc(
    f::algebraic, 
    X0::(name=realcons),
    { 
        digits::posint:= Digits,
        tolerance::positive:= .5*10^(2-Digits),
        maxiters::posint:= 20,
        criterion::identical(relative, absolute, residual):= ':-relative'
    }
)
local 
    X:= Array(1..2, [rhs(X0)+1, rhs(X0)], datatype= sfloat),
    x:= lhs(X0),
    F:= subs(
        _f= unapply(f,x), 
        proc(x) option remember; evalf(_f(x)) end proc
    ),
    N:= subs(_d= unapply(diff(f,x), x), x-> x-F(x)/evalf(_d(x))),
    Crit:= table([
        ':-absolute'= ((a,b,t)-> abs(a-b) < t),
        ':-relative'= ((a,b,t)-> abs(1-b/a) < t),
        ':-residual'= ((a,b,t)-> abs(F(a)) < t)
    ]),
    k
;
    Digits:= digits;
    for k from 3 to 2+maxiters do
        if Crit[criterion](X[k-1], X[k-2], tolerance) then break fi;
        X(k):= N(X[k-1])
    od;
    if k > 2+maxiters then print("Did not converge") fi;
    seq(X[2..])
end proc
:   
f:= sin(2*x)+3*x+2*cos(x)^2-61:
Newton(f, x=0, digits= 15, tolerance= .5e-13, criterion= relative);

@vsnyder I meant a transcript of the Fortran calls to EvalMapleExpression, not a transcript of how the results appear  within Maple. But at this point I think that it's clear to both of us what you were doing wrong: The procedure definition, or any expression, must appear as a single call to EvalMapleExpression. It's make no difference whether the argument to that call contains line breaks. Analogously, if you were entering code into Maple by hand, regardless of interface, no single expression (even a thousand-line module) could be split over multiple execution groups (command prompts); nor would you be allowed to press Return until it was completely entered. 

The file attachment tool in MaplePrimes is the green uparrow rather than the more-usual paperclip icon on a lot of sites. At least the bright color makes it easier to see.

@Rookieplayer Yes, I noticed many SLATEC routines related to piecewise polynomials (just use your browser's in-page search feature (usually Ctrl-F) on the page https://gams.nist.gov/cgi-bin/serve.cgi/PackageModules/SLATEC), but nothing for non-polynomials. Of course, your functions are piecewise-analytic, so they could be approximated by piecewise polynomials; however, I think that the approach that I described in my previous Reply will be far more accurate, far more computationally efficient, and much easier to implement.

By the way, in case you don't already know this: You should check out Maple's command CodeGeneration:-Fortran, which I used to generate HVS above.

@Rookieplayer I can't find piecewise in SLATEC either. Before sending your code to Fortran, convert all the piecewises to Heaviside (this could be done all at once with the subsindets command). I can't find Heaviside in SLATEC either, but it'd be trivial to implement. Substitute HVS (or another Fortran-acceptable name) for Heaviside (this could be done all at once with subs). Then use Fortran code such as this for HVS:

doubleprecision function HVS (x)
    doubleprecision x
    if (x .le. 0.0D0) then
         HVS = 0.0D0
         return
    else
         HVS = 0.1D1
         return
    end if
end

My latex(...) results are quite different from yours:

restart:
Physics:-Version();
         The "Physics Updates" package is not installed

interface(version);
Standard Worksheet Interface, Maple 2021.0, Windows 10, March 5 
   2021 Build ID 1523359

interface(typesetting);
                            extended
interface(prettyprint);
                               2
latex:-Settings();
[cacheresults = true, commabetweentensorindices = false, 
  invisibletimes = " ", leavespaceafterfunctionname = false, 
  linelength = 66, powersoftrigonometricfunctions = mixed, 
  spaceaftersqrt = true, usecolor = true, 
  usedisplaystyleinput = true, useimaginaryunit = I, 
  useinputlineprompt = true, userestrictedtypesetting = false, 
  usespecialfunctionrules = true, 
  usetypesettingcurrentsettings = false]

sol:= u(r,t) = 
    invlaplace(
        BesselJ(0,10*(-s)^(1/2)*r)/BesselJ(0,20*(-s)^(1/2))*s/(s^2+1),
        s, t
     ) 
     - invlaplace(
           BesselJ(0,10*(-s)^(1/2)*r)/BesselJ(0,20*(-s)^(1/2))/s,
           s, t
     )
     - cos(t) + 1
:
latex(sol);
u \! \left(r , t\right) = 
\mathcal{L}^{-1}\! \left(\frac{J_{0}\! \left(10 \sqrt{-s}\, r \right)
 s}{J_{0}\! \left(20 \sqrt{-s}\right) \left(s^{2}+1\right)}, s , t\right)
-\mathcal{L}^{-1}\! \left(\frac{J_{0}\! \left(10 \sqrt{-s}\, r \right)}
{J_{0}\! \left(20 \sqrt{-s}\right) s}, s , t\right)-\cos \! \left(t \right)+1