Hi,

I need to run the following procedure a couple of million times. Although it works, Maple sometimes chokes for no apparent reason (if there is a reason, please let me know). I was wondering whether an expert could help me tweak the procedure (or possibly rewrite it) to achieve the best possible performance. I am planning to use Grid:-Map or, if possible, Threads:-Map.

 

generateNonlinearModelsPlus := proc(model::list,fullmodel::list,vars::list:=[x,y,z])
description "This function generates a list of all models with one more monomial from the full model":
local tab::table(),n:=nops(model),i,j,k:=1,ans,terms,aaa,allmoncoefThreads:
# local procedure
allmoncoefThreads := proc(f::list,vars::list)
description "This function finds the monomials multipled by their coefficients for each expression (equation) of a list.":
local n:=numelems(f),i,mon:=[seq](0,i=1..n),M,cc:=[seq](0,i=1..n),ans:
for i from 1 to n do
  cc[i]:=[coeffs](expand(f[i]),vars, 'M'):
  mon[i]:=[M]:
end do:
ans:=[seq](zip((ww,vv)->ww*vv,cc[i],mon[i]),i=1..n):
return(ans)
end proc:
# main part
ans:=zip((w,v)->expand(simplify(v-w)),model,fullmodel): # Find the monomials that are not in model
terms:=allmoncoefThreads(ans,vars): # Separate the monomials
#
for i from 1 to n do
   aaa:=model:
   for j from 1 to nops(terms[i]) do
       aaa[i]:=model[i]+terms[i,j]:
       tab[k]:=aaa:
       k:=k+1:
   end do:
end do:
tab:=convert(tab,list):
return(tab):
end proc:

Here is an example of how I run it: 

model:=[y*alpha[1, 2], z*alpha[2, 3], x^3*alpha[3, 10] + x*alpha[3, 1] + alpha[3, 0]]:

fullmodel:=[x^3*alpha[1, 10] + x^2*y*alpha[1, 11] + x^2*z*alpha[1, 12] + x*y^2*alpha[1, 13] + x*y*z*alpha[1, 14] + x*z^2*alpha[1, 15] + y^3*alpha[1, 16] + y^2*z*alpha[1, 17] + y*z^2*alpha[1, 18] + z^3*alpha[1, 19] + x^2*alpha[1, 4] + x*y*alpha[1, 5] + x*z*alpha[1, 6] + y^2*alpha[1, 7] + y*z*alpha[1, 8] + z^2*alpha[1, 9] + x*alpha[1, 1] + y*alpha[1, 2] + z*alpha[1, 3] + alpha[1, 0], x^3*alpha[2, 10] + x^2*y*alpha[2, 11] + x^2*z*alpha[2, 12] + x*y^2*alpha[2, 13] + x*y*z*alpha[2, 14] + x*z^2*alpha[2, 15] + y^3*alpha[2, 16] + y^2*z*alpha[2, 17] + y*z^2*alpha[2, 18] + z^3*alpha[2, 19] + x^2*alpha[2, 4] + x*y*alpha[2, 5] + x*z*alpha[2, 6] + y^2*alpha[2, 7] + y*z*alpha[2, 8] + z^2*alpha[2, 9] + x*alpha[2, 1] + y*alpha[2, 2] + z*alpha[2, 3] + alpha[2, 0], x^3*alpha[3, 10] + x^2*y*alpha[3, 11] + x^2*z*alpha[3, 12] + x*y^2*alpha[3, 13] + x*y*z*alpha[3, 14] + x*z^2*alpha[3, 15] + y^3*alpha[3, 16] + y^2*z*alpha[3, 17] + y*z^2*alpha[3, 18] + z^3*alpha[3, 19] + x^2*alpha[3, 4] + x*y*alpha[3, 5] + x*z*alpha[3, 6] + y^2*alpha[3, 7] + y*z*alpha[3, 8] + z^2*alpha[3, 9] + x*alpha[3, 1] + y*alpha[3, 2] + z*alpha[3, 3] + alpha[3, 0]]:

vars:=[x,y,z]:

ans:=generateNonlinearModelsPlus(model,fullmodel,vars)

Many thanks.

For quite some time, I have wanted to solve the system attached in "test" using Maple. The smallest solution in natural numbers x, y, and z test.mw

>  (1) >  (2) (3) (4)

Download test.mw

is known, and all these numbers are less than 4 × 10¹². Is this possible in Maple?

(x=1633780814400; y=252782198228; z=3474741058973)

M being some sparse matrix, plots:-sparsematrixplot(M) displays a patchwork of squares with two different colors: one for elements of M different 0 and the other for elements of M equal to 0.
The squares are separated by what visually appears like blank lines (in fact tere are no such lines, just gaps between the different squares making them not adjacent).

So my question: Is it possible, without resorting to any workaround based upon PLOT(POLYGONS(..)) or plottools functions (I know how to write such workarounds and I'm not interested in that) to suppress these "white lines", or otherwise said to make the squares adjacent? 
I thought that maybe some (hidden?) option acting like style=patchnogrid could exist, or to an undocumented feaure...

Thanks for your reply.

The title is the question. I am fairly certain that Maple in itself (up to 2023) does not have this. There is of course Statistics:-Histogram, but that is a 1-D histogram only.

Given 2 Vectors of floats, say, xvec[n] and yvec[n]. The idea would be to divide up a 2-D grid into "pixels" that accumulate a count when (xvec[i],yvec[i]) falls into the range for a pixel (xrange,yrange). Plotting each of these pixels as rectangular column creates the histogram.

I could of course program one myself, but as I am lazy and this seems like a common kind of plot (at least in statistics and in physics), someone may well have done such a thing. DDG search did not find anything, at least not wrt. Maple.

Here is an example plot close to what I am thinking of (random plot pulled off the web). The color scheme of the plot is not what I am after (although it does look nice).

Here is a snippet of code that returns an error message. The details of the h procedure should be irrelevant.

h:=proc(k,x)
local z,w;
nans(parameters=[k,x]);
z:=nans(0.98);
w:=rhs(z[2]),rhs(z[3])
end proc;
 h := proc (k, x) local z, w; nans(parameters = [k, x]); z := 

    nans(.98); w := rhs(z[2]), rhs(z[3]) end proc;

h(7.7,0.3);
0.310903619302454933005944397083, 0.731944275983583692659532399302

plot([h(7.7,x),x=0..30]);
Error, (in dsolve/numeric/process_parameters) 'parameters' must be specified as a list of numeric values

Note that the line h(7.7,0.3); is just there to show that h functions well. It accepts 2 numerical inputs and spits out a pair of numerical values.

Why is plot not plotting? Why is plot asking about inner workings of h?

If I try

f:=proc(x,y)

cos(x),sin(x*y) 
end proc;
plot([f(x,7.7),x=0..3]);
that woks fine. Plot does not inquire about f.

I am trying to select from a default list and a modified list to get a third combined list of the unmodified default elements and the mdfified elements plus extra elements. This is for plot data. I have managed to extract the colour data using something I did a couple of years ago. Though at this stage I don't relly know how that works either but it works.

>  restart >  with(ListTools) (1) >  PlotDefaults:=['colour' = ':-blue', symbol = ':-solidcircle', ':-symbolsize' = 8,thickness=2]; #for points and lines (2) >  Inputs:=['color'=[red, black, blue],symbol=square,':-linestyle'=dash]; #'color'=[red, black, blue], >    (3) >  if has(Inputs,{colour,color}) then Colourlist:=remove(has,remove(has,Flatten(eval([':-color',':-colour'],Inputs)),':-colour'),':-color') else Colourlist:=remove(has,remove(has,Flatten(eval([':-color',':-colour'],PlotDefaults)),':-colour'),':-color'); end if (4) Make this list >  plotdata:=[symbol=square,':-linestyle'=dash,':-symbolsize' = 8,thickness=2] (5) >  mx:=max(nops(PlotDefaults),nops(Inputs)) (6) >  Plotdata:=[]; for i to mx do end do >    >    (7)

Download 2026-05-18_Q_Select_from_Two_Lists_to_get_New_List.mw

Currently I have Maple versions 2023,2025, and 2026 installed on Windows 11. Today I installed a workbook package containing a module that I just completed using the PackageTools installer in Maple 2026.. To my surprise, I found that a package installed from Maple 2026 was also available in Maple 2023 and, conversely, a package installed in Maple 2023 was automatically available in Maple 2026. i noticed that, with the exception of the Maple Customer Support Updates, the toolbox directory is no longer broken down by versions. I also noticed that the directory containing the module installed by Maple 2026 was named by the workbook instead of the module name (ie. hopfwords.maple). As I recall, the toolboxes used to be version dependent. 

The question is to what extent can one assume that a package created in Maple 2026 will be compatible with at least the more recent versions of Maple, I am also wondering why the directory name is now the workbook name instead of the module name. 

restart;

with(plots): with(LinearAlgebra):

 

# TFSB Coefficients (symbolic in u)

beta0 := u -> (sin(u)*u^3 - 12*u^2 - 24*cos(u) + 24)/(12*(sin(u)*u + 2*cos(u) - 2)*u^2):

beta1 := u -> (5*sin(u)*u^3 + 12*cos(u)*u^2 + 24*cos(u) - 24)/(6*(sin(u)*u + 2*cos(u) - 2)*u^2):

beta2 := u -> beta0(u):

rho0 := u -> ((-u^2-12)*cos(u) - 5*u^2 + 12)/(12*(sin(u)*u + 2*cos(u) - 2)*u^2):

rho1 := u -> (-7*cos(u)*u^3 + 27*sin(u)*u^2 + 120*sin(u) - 120*u)/(60*u^2*(cos(u)*u + 2*u - 3*sin(u))):

rho2 := u -> -rho0(u):

 

# Secondary coefficients (simplified versions)

beta00 := u -> 13/42 - 9*u^2/7840:

beta10 := u -> 1/6 + u^2/720:

beta20 := u -> 1/42 - 17*u^2/70560:

beta01 := u -> 187/1680 + 611*u^2/705600:

beta11 := u -> 11/30 - 29*u^2/25200:

beta21 := u -> 37/1680 + 67*u^2/235200:

beta02 := u -> 11/70 + 491*u^2/352800:

beta12 := u -> 9/10 - 31*u^2/8400:

beta22 := u -> 31/70 + 811*u^2/352800:

 

rho01 := u -> 2/105 + 407*u^2/1058400:

rho11 := u -> -19/210 + 41*u^2/105840:

rho21 := u -> -1/168 - 101*u^2/529200:

rho02 := u -> 53/1680 + 1633*u^2/2116800:

rho12 := u -> 8/105 - 4*u^2/6615:

rho22 := u -> -101/1680 - 2273*u^2/2116800:

 

# Problem definition

omega := 1:

epsilon := 3*Pi/2:

phi := x -> 3*sin(x) - 5*cos(x):  # history function

 

f := (x, v, vp, vd) -> -v - vd + 3*cos(x) + 5*sin(x):

g := proc(x, v, vp, vd, vdp)

    local fx, fv, fvp, fvd;

    fx := -3*sin(x) + 5*cos(x);

    fv := -1;

    fvp := 0;

    fvd := -1;

    return fx + fv*vp + fvp*0 + fvd*vdp;

end proc:

 

# Initial conditions

a := 0: b := 10:

v0 := -5: vp0 := 3:

 

# Variable step-size parameters

tol := 1e-10:

h_min := 0.01:

h_max := 0.5:

h_init := Pi/8:

 

# Store results

X := [a]: V := [v0]: Vp := [vp0]:

h_curr := h_init:

x_curr := a:

v_curr := v0:

vp_curr := vp0:

 

# For history: need v at x-epsilon

get_v_delayed := proc(xx)

    if xx < a then return phi(xx);

    else

        # Interpolate from stored solution

        idx := 1;

        while idx < nops(X) and X[idx] < xx do idx := idx+1; end do;

        if idx = 1 then return phi(xx);

        elif X[idx] = xx then return V[idx];

        else

            # Linear interpolation

            return V[idx-1] + (V[idx]-V[idx-1])*(xx-X[idx-1])/(X[idx]-X[idx-1]);

        end if;

    end if;

end proc:

 

# Newton solver for block

solve_block := proc(x0, v0, vp0, h, omega)

    local u, bet0, bet1, bet2, rho0, rho1, rho2, bet00, bet10, bet20, bet01, bet11, bet21, bet02, bet12, bet22,

          rho01, rho11, rho21, rho02, rho12, rho22, F, J, V0, V1, V2, Vp0, Vp1, Vp2, tolN, iter, dv, dV;

   

    u := omega*h;

    bet0 := beta0(u); bet1 := beta1(u); bet2 := beta2(u);

    rho0 := rho0(u); rho1 := rho1(u); rho2 := rho2(u);

    bet00 := beta00(u); bet10 := beta10(u); bet20 := beta20(u);

    bet01 := beta01(u); bet11 := beta11(u); bet21 := beta21(u);

    bet02 := beta02(u); bet12 := beta12(u); bet22 := beta22(u);

    rho01 := rho01(u); rho11 := rho11(u); rho21 := rho21(u);

    rho02 := rho02(u); rho12 := rho12(u); rho22 := rho22(u);

   

    # Initial guesses

    V1 := v0 + h*vp0;

    V2 := v0 + 2*h*vp0;

    Vp1 := vp0;

    Vp2 := vp0;

   

    tolN := 1e-12;

    for iter from 1 to 10 do

        # Compute delayed values

        vd0 := get_v_delayed(x0 - epsilon);

        vd1 := get_v_delayed(x0 + h - epsilon);

        vd2 := get_v_delayed(x0 + 2*h - epsilon);

        vdp0 := (get_v_delayed(x0 - epsilon + 1e-8) - vd0)/1e-8;

        vdp1 := (get_v_delayed(x0 + h - epsilon + 1e-8) - vd1)/1e-8;

        vdp2 := (get_v_delayed(x0 + 2*h - epsilon + 1e-8) - vd2)/1e-8;

       

        # Compute gamma and g

        gam0 := f(x0, v0, vp0, vd0);

        gam1 := f(x0+h, V1, Vp1, vd1);

        gam2 := f(x0+2*h, V2, Vp2, vd2);

        g0 := g(x0, v0, vp0, vd0, vdp0);

        g1 := g(x0+h, V1, Vp1, vd1, vdp1);

        g2 := g(x0+2*h, V2, Vp2, vd2, vdp2);

       

        # Residuals

        F1 := h*vp0 - (V1 - v0 + h^2*(bet00*gam0 + bet10*gam1 + bet20*gam2)

              + h^3*(rho01*g0 + rho11*g1 + rho21*g2));

        F2 := h*Vp1 - (V1 - v0 + h^2*(bet01*gam0 + bet11*gam1 + bet21*gam2)

              + h^3*(rho01*g0 + rho11*g1 + rho21*g2));

        F3 := h*Vp2 - (V1 - v0 + h^2*(bet02*gam0 + bet12*gam1 + bet22*gam2)

              + h^3*(rho02*g0 + rho12*g1 + rho22*g2));

        F4 := V2 - (2*V1 - v0 + h^2*(bet0*gam0 + bet1*gam1 + bet2*gam2)

              + h^3*(rho0*g0 + rho1*g1 + rho2*g2));

       

        F := Vector([F1, F2, F3, F4]);

        if LinearAlgebra:-Norm(F) < tolN then break; end if;

       

        # Approximate Jacobian (finite differences)

        J := Matrix(4,4);

        delta := 1e-6;

        for j from 1 to 4 do

            V_pert := Vector([V1, V2, Vp1, Vp2]);

            V_pert[j] := V_pert[j] + delta;

            V1p := V_pert[1]; V2p := V_pert[2]; Vp1p := V_pert[3]; Vp2p := V_pert[4];

            gam1p := f(x0+h, V1p, Vp1p, get_v_delayed(x0+h-epsilon));

            gam2p := f(x0+2*h, V2p, Vp2p, get_v_delayed(x0+2*h-epsilon));

            g1p := g(x0+h, V1p, Vp1p, get_v_delayed(x0+h-epsilon),

                     (get_v_delayed(x0+h-epsilon+1e-8)-get_v_delayed(x0+h-epsilon))/1e-8);

            g2p := g(x0+2*h, V2p, Vp2p, get_v_delayed(x0+2*h-epsilon),

                     (get_v_delayed(x0+2*h-epsilon+1e-8)-get_v_delayed(x0+2*h-epsilon))/1e-8);

           

            F1p := h*vp0 - (V1p - v0 + h^2*(bet00*gam0 + bet10*gam1p + bet20*gam2p)

                   + h^3*(rho01*g0 + rho11*g1p + rho21*g2p));

            F2p := h*Vp1p - (V1p - v0 + h^2*(bet01*gam0 + bet11*gam1p + bet21*gam2p)

                   + h^3*(rho01*g0 + rho11*g1p + rho21*g2p));

            F3p := h*Vp2p - (V1p - v0 + h^2*(bet02*gam0 + bet12*gam1p + bet22*gam2p)

                   + h^3*(rho02*g0 + rho12*g1p + rho22*g2p));

            F4p := V2p - (2*V1p - v0 + h^2*(bet0*gam0 + bet1*gam1p + bet2*gam2p)

                   + h^3*(rho0*g0 + rho1*g1p + rho2*g2p));

           

            Fp := Vector([F1p, F2p, F3p, F4p]);

            J[1..4, j] := (Fp - F)/delta;

        end do;

       

        dV := LinearAlgebra:-LinearSolve(J, -F);

        V1 := V1 + dV[1]; V2 := V2 + dV[2]; Vp1 := Vp1 + dV[3]; Vp2 := Vp2 + dV[4];

    end do;

   

    return [V1, V2, Vp1, Vp2];

end proc:

 

# Main variable step-size loop

printf("Variable step-size integration for Example 1\n");

printf("tol = %e, h_init = %f\n", tol, h_init);

 

while x_curr < b - 1e-12 do

    # Try current step

    sol := solve_block(x_curr, v_curr, vp_curr, h_curr, omega);

    V1 := sol[1]; V2 := sol[2]; Vp1 := sol[3]; Vp2 := sol[4];

   

    # Compute with two half-steps

    sol_half1 := solve_block(x_curr, v_curr, vp_curr, h_curr/2, omega);

    V_mid := sol_half1[2]; Vp_mid := sol_half1[4];

    sol_half2 := solve_block(x_curr + h_curr/2, V_mid, Vp_mid, h_curr/2, omega);

    V2_half := sol_half2[2];

   

    # Error estimate

    err := abs(V2 - V2_half) / (2^6 - 1);

   

    if err < tol then

        # Accept step

        x_next := x_curr + 2*h_curr;

        X := [op(X), x_curr + h_curr, x_next];

        V := [op(V), V1, V2];

        Vp := [op(Vp), Vp1, Vp2];

        x_curr := x_next;

        v_curr := V2;

        vp_curr := Vp2;

        printf("x = %7.4f, h = %8.5f, err = %12.5e\n", x_curr, h_curr, err);

       

        # Adjust step size

        if err < tol/2 then

            h_curr := min(2*h_curr, h_max);

        end if;

    else

        # Reject step, reduce h

        h_curr := max(h_curr/2, h_min);

        printf("  Rejecting, new h = %8.5f\n", h_curr);

    end if;

end do:

 

# Exact solution for comparison

exact := x -> 3*sin(x) - 5*cos(x);

errors := [seq(abs(V[i] - exact(X[i])), i=1..nops(X))];

 

# Visualization

p1 := pointplot([seq([X[i], errors[i]], i=1..nops(X))], color=red, symbol=circle,

                title="Example 1: Variable Step-Size TFSB - Absolute Errors",

                labels=["x", "Error"], labeldirections=[horizontal,vertical]);

p2 := plot([[x_curr, h_curr]], x=a..b, style=point, color=blue,

            title="Step-size evolution", labels=["x", "h"]);

display(p1);

display(p2);

 

printf("\nFinal results for Example 1:\n");

printf("Number of steps: %d\n", nops(X)-1);

printf("Maximum error: %e\n", max(errors));

printf("Final step-size: %f\n", h_curr);

Since Maple version 2026, I have noticed that the orientation option with the contourplot3d command has no effect.

It is easy to control by using the example provided in the online help and adding, for example, orientation=[20,10,10].

Thank you for your help.

Best regards.

...the essence of plane geometry is hidden within the following puzzle:
Given is a closed curve C. It is assumed to be non-self-intersecting, convex, and continuously differentiable everywhere (a closed Jordan curve). Let line segment AB be a chord of this curve, having a fixed length l. A point P lies on this chord at a fixed distance a from A and b from B, such that l = a + b. An orientation (or direction of circulation) is now assigned to the curve. The chord is then moved continuously along the closed curve in this assigned direction of circulation. As it moves, point P traces out a so-called locus curve O, which—upon completion of one full revolution of the chord—also forms a closed curve lying entirely within C.
The task is to calculate the area of ​​the region between C and O (i.e., the area lying inside C but outside O). Divide the result by the product of a and b, and then apply the "identify" function to the outcome.

On my system (Windows 11) I can only remove an entry from the favorites when a document is open. With all documents closed removing does not work.

Can someone confirm?

I was using pdsolve to solve for a 1-D longitudinal wave equation.  My particular problem added an external stimulus that is not mechanical in original, so the acoustic wave velocity has to be modeled as a variable c(H) and I also have to add a du/dx term in the model.  You can see my question posted on April 17, 2026.  I was told that pdsolve does not handle such a problem.  The suggestion in the posting is to use a finite difference method.  

I am verifying that approach by solving a simplier problem where c is constant but with an initial uniformly stretched material.  The solution does not seem to be physical.  The material should relax throughout the whole length of the material, but the solution shows relaxation at the end and stay uniformly stretched at the center.  I ran the problem with pdsolve and get a different result that I think is more realistic.

Is there something I can tweak in the finite difference approach to overcome that issue?

pde_finite_difference_method_linear_ic.mw

pdsolve_exercise_damping_ini_linear_a.mw

Is there a workaround to being able to call object method inside the debugger? I use Objects everywhere and hard to debug things when not able to try calls inside the debugger.

Here is an example.

This worksheet below has the code to reproduce this. Ofcourse the code works when I run it in the debugger using "continue" button, or using the "next" or "step". But I need sometimes to check something before doing this.

>  interface(version); >  restart; >  kernelopts('assertlevel'=2): >  module person_type()     option object;     local _name::string:="";         export ModuleCopy::static := proc( _self::person_type, proto::person_type, name::string)       _name:=name;       NULL;     end proc;     export get_name::static:=proc(_self)        _self:-_name;     end proc; end module; >  foo:=proc(o::person_type)::string;   local name;   DEBUG();   name:=o:-get_name();   name; end proc; >  o:=Object(person_type,"me"); >  o:-get_name(); >  #run this. Debugger will now open foo(o);

Download why_calling_object_fail_in_debugger.mw

Please see the attached worksheet:

Question 1:  How to plot a partial portion of an array s2 using "dataplot," say the first 50 elements.

Question 2:  I have also tried to use "plot" to plot the first 50 elements s2 and keep on getting errors or strange looking plot with horizontal lines.  I have tried plot(s2,1..50), plot(s2,x=1..50), plot(s2,m=1..50), plot(s2[m],m=1..50), plot(s2[1..50],1..50], plot(1..50,s2[1..50]).

plot_exercise_a.mw

How i can apply this transfromation which is special and i can't undrestand how from that solution he did thus changes  and get the solution of ode  ,  i reached untill eq(3.5) but after that i can't  figure out what he did   and how we get eq(3.7)-(3.9) (3.12),(3.13) and how we can   tranform solution  for the  eq(3.5)

fp.mw