MaplePrimes Questions

What does Error, (in dsolve/numeric/bvp) bad index into Matrix mean?
Also, I'm trying to run it, it is slow, any suggestions?

restart;
with(Student[VectorCalculus]);
with(DynamicSystems);
with(DEtools);
with(PDEtools, ReducedForm, declare, diff_table, dsubs);
NULL;
 #Digits:= trunc(evalhf(Digits)); #generally a very efficient setting

# Here we solve a 1D problem in 3 regions. In each region, we have concentration and potential (c,phi) distributions,
# We first solve the unperturbed steady-state problem and then the linearized perturbation problem (which rely on the steady state).
# Each region is defined in x = 0..1, and the regions are connected by interface conditions that 
# connect (c1(1),phi1(1)) to (c2(0),phi2(0)) and (c2(1),phi2(1)) to (c3(0),phi3(0))

Q:=10;   omega:=100;     J0:= 0.01;   # parameters
                            Q := 10

                          omega := 100

                           J0 := 0.01

# The unperturbed steady-state

c1:=1-J0/2*x: 
c3:=1-J0/2*(x-1):  
c12:= eval(c1,x=1); 
c32 := eval(c3,x=0); 
S1:=sqrt(Q^2+4*c12^2): 
S3:=sqrt(Q^2+4*c32^2):  
c21:=eval((S1-Q)/2); 
c23:=eval((S3-Q)/2);  
I0:=fsolve(Q*i0/2/J0*ln((J0*S1-Q*i0)/(J0*S3-Q*i0))=(J0-S1+S3)/2,i0);  
V:=(I0/J0+1)*ln(c32/c12)+ln((c21+Q)/(c23+Q))+(J0+2*c23-2*c21)/Q; # the potential drop across the system 
c2:=solve(y-c21+Q*I0/2/J0*ln((Q*I0-Q*J0-2*J0*y)/(Q*I0-Q*J0-2*J0*c21))=-J0/2*x,y):  
phi1:=I0/J0*ln(c1)+V: 
phi3:=I0/J0*ln(c3): 
dphi1:=diff(phi1,x); 
dphi3:=diff(phi3,x); 
phi21:=I0/J0*ln(c12)+V-0.5*ln((c21+Q)/c21); 
phi2:=(2*c21-2*c2+Q*phi21-J0*x)/Q: 
dphi2:=diff(phi2,x); 
dphi12 := eval(dphi1, x=1); 
dphi21 := eval(dphi2, x=0); 
dphi23 := eval(dphi2, x=1); 
dphi32 := eval(dphi3, x=0); 
INT1 := int(1/(2*c1), x = 0 .. 1); 
INT2 := int(1/(2*c2 + Q), x = 0 .. 1); 
INT3 := int(1/(2*c3), x = 0 .. 1); 
INT := INT1 + INT2 + INT3;
                      c12 := 0.9950000000

                       c32 := 1.005000000

                      c21 := 0.09804129000

                      c23 := 0.1000024500

                      I0 := 0.01419804328

                       V := 0.02539628566

                              0.007099021640   
                dphi1 := - --------------------
                           1 - 0.005000000000 x

                              0.007099021640        
           dphi3 := - ------------------------------
                      1.005000000 - 0.005000000000 x

                     phi21 := -2.299074561

dphi2 := (0.001000000000 LambertW(-0.2818670588 exp(-0.2818670588

   - 0.0007043224058 x)))/(1

   + LambertW(-0.2818670588 exp(-0.2818670588 - 0.0007043224058 x)

  )) - 0.001000000000


                   dphi12 := -0.007134695118

                   dphi21 := -0.001392499832

                   dphi23 := -0.001391964358

                   dphi32 := -0.007063703124

                      INT1 := 0.5012541824

                     INT2 := 0.09805801917

                      INT3 := 0.4987541511

                       INT := 1.098066353


sys1 := {
-omega*C11(x)+diff(diff(C12(x), x), x)=0,
omega*C12(x)+diff(diff(C11(x), x), x) = 0,
-omega*C21(x)+diff(diff(C22(x), x)+(c2*sigma2-C22(x)*dphi2*Q)/(2*c2+Q), x) =0,
 omega*C22(x)+diff(diff(C21(x), x)+(c2*sigma1-C21(x)*dphi2*Q)/(2*c2+Q), x) = 0,
-omega*C31(x)+diff(diff(C32(x), x), x)=0,
omega*C32(x)+diff(diff(C31(x), x), x) = 0
}:

sys2 := {
diff(FA1(x), x) = C11(x)*dphi1/c1,
diff(FA2(x), x) = C21(x)*dphi2/(c2+Q/2),
diff(FA3(x), x) = C31(x)*dphi3/c3,
diff(FB1(x), x) = C12(x)*dphi1/c1,
diff(FB2(x), x) = C22(x)*dphi2/(c2+Q/2),
diff(FB3(x), x) = C32(x)*dphi3/c3
}: 

Bc := {
C11(0) = 0, C12(0) = 0,  C31(1) = 0, C32(1) = 0,
FA1(0) = 0, FB1(0) = 0,  FA3(1) = 0, FB3(1) = 0, 

2*C11(1)/c12 = C21(0)/(c21+Q)+C21(0)/c21, 
2*C12(1)/c12 = C22(0)/(c21+Q)+C22(0)/c21,
C21(1)/(c23+Q)+C21(1)/c23 = 2*C31(0)/c32,
C22(1)/(c23+Q)+C22(1)/c23 = 2*C32(0)/c32,

D(C11)(1)+dphi12*C11(1)-sigma1/2-c12*D(FA1)(1) = D(C21)(0)+dphi21*C21(0)-(c21+Q)*sigma1/(2*c21+Q)-(c21+Q)*D(FA2)(0),
D(C12)(1)+dphi12*C12(1)-sigma2/2-c12*D(FB1)(1) = D(C22)(0)+dphi21*C22(0)-(c21+Q)*sigma2/(2*c21+Q)-(c21+Q)*D(FB2)(0),
D(C11)(1)-dphi12*C11(1)+sigma1/2+c12*D(FA1)(1) = D(C21)(0)-dphi21*C21(0)+c21*sigma1/(2*c21+Q)+c21*D(FA2)(0),
D(C12)(1)-dphi12*C12(1)+sigma2/2+c12*D(FB1)(1) = D(C22)(0)-dphi21*C22(0)+c21*sigma2/(2*c21+Q)+c21*D(FB2)(0),

D(C31)(0)+dphi32*C31(0)-sigma1/2-c32*D(FA3)(0) = D(C21)(1)+dphi23*C21(1)-(c23+Q)*sigma1/(2*c23+Q)-(c23+Q)*D(FA2)(1),
D(C32)(0)+dphi32*C32(0)-sigma2/2-c32*D(FB3)(0) = D(C22)(1)+dphi23*C22(1)-(c23+Q)*sigma2/(2*c23+Q)-(c23+Q)*D(FB2)(1),
D(C31)(0)-dphi32*C31(0)+sigma1/2+c32*D(FA3)(0) = D(C21)(1)-dphi23*C21(1)+c23*sigma1/(2*c23+Q)+c23*D(FA2)(1),
D(C32)(0)-dphi32*C32(0)+sigma2/2+c32*D(FB3)(0) = D(C22)(1)-dphi23*C22(1)+c23*sigma2/(2*c23+Q)+c23*D(FB2)(1)
}:
 
 


all_sys := sys1 union sys2 union Bc:
sol1 := dsolve(all_sys, initmesh = 100, maxmesh = 15000, numeric, method = bvp[midrich], output = listprocedure):
#(all_sys, numeric, method = bvp[midrich]);

Error, (in dsolve/numeric/bvp) bad index into Matrix

Hello,

 I would like to bulild a Gaussian Mixture Model in Maple using the most straightforward apprach:

Apperantly something goes wrong, PDF in (5) is a product not a sum and variance is quadratic in w (6) not linear as it should be.

How to correct this?

Helllo folkx :)

What is wrong in my functions f and g ?

I don't catch the point.

Thanks.

JM--

interface(version);
    Standard Worksheet Interface, Maple 2025.2, Windows 11,
     November 11 2025 Build ID 1971053

with(RandomTools):

Generate(float(range = 2.532000 .. 7.723000, digits = 4));
                             2.537

f := x -> Generate(float(range = 0 .. 100)):

f(30);
                            94.3141

f(0.200000);
                            23.9622

randomize();
                         2413841372727

f(30);
                            86.2657

plot(f(x), x = 1 .. 100.000000);

plots[listplot](sort([seq(Generate(float(range = 0.032100 .. 162.000000, digits = 3)), i = 1 .. 20)]));

plots[listplot](sort([seq(Generate(float(range = 0.032100 .. 162.000000, digits = 3, method = uniform)), i = 1 .. 20)]));

g := ((x -> Generate(float(range = 0 .. x))) assuming (x*float and 0 < x)):

plot(g(x), x = 1 .. 50);
Error, (in RandomTools:-Generate) invalid input: unknown expects value for keyword parameter range to be of type numeric .. numeric, but received 0 .. x

Is the following behavior expected? using  t=0 .. 3*Pi  vs.  t=0 .. round(3*Pi) give different plots.  Why? Should not the plot be the same?

ps. do not know if I asked this before or not. I can't remember. Search in Mapleprimes is not easy.

interface(version);

`Standard Worksheet Interface, Maple 2025.2, Windows 10, November 11 2025 Build ID 1971053`

restart;

ode := diff(y(t), t$2) + y(t)=0;
DEtools:-DEplot(ode, y(t), t=0 .. 3*Pi, y=-1 .. 1,[[y(0)=1,D(y)(0)=0]],linecolor=blue);

diff(diff(y(t), t), t)+y(t) = 0

DEtools:-DEplot(ode, y(t), t=0 .. round(3*Pi), y=-1 .. 1,[[y(0)=1,D(y)(0)=0]],linecolor=blue);

 

 

Download deplot_with_round_jan_17_2026.mw

I am looking for ways to clear remember tables in a Maple session to ensure that a Maple command/procedure is executed as if it were being called for the first time. I currently use

map(forget, [anames()])

but I am not sure if that statement clears all existing remeber tables (i.e. including system remember tables)?

Are there other ways to do this? Restart is not an option since it clears almost everything.

Edit: This question is not about disabling remember tables.

I'm looking for the general solution to the attached differential equation. Maple doesn't provide it. What am I doing wrong?

restart

ode5 := diff(y(x), x) = (8*y(x)*b-32*b^2*x/y(x)-64*b^2*x^2*y(x))/(3*y(x)^2+8*b*x-16*b^2*x^2/y(x)^2)

diff(y(x), x) = (8*y(x)*b-32*b^2*x/y(x)-64*b^2*x^2*y(x))/(3*y(x)^2+8*b*x-16*b^2*x^2/y(x)^2)

(1)

simplify(ode5)

diff(y(x), x) = ((64*b^2*x^2-8*b)*y(x)^3+32*b^2*x*y(x))/(-3*y(x)^4-8*b*x*y(x)^2+16*b^2*x^2)

(2)

dsolve(ode5, y(x))

NULLNULL

Download testdgl5.mw

Hi,

I am using the following (dummy) code. I would appreciate any help generating a plot with multiple functions and a legend. I am getting an error message stating that the parameter "v" is unknown. I believe the issue is due to the complexity of the functions (fun1 and fun2) that use the "piecewise" command in my code.  

fun1 := piecewise(cond1, a(x,v), cond2, b(x,v))

fun2 := piecewise(cond1, c(x,v), cond2, d(x,v))

with(plots):
myPlotFunction := proc(v)
local p1, p2:
p1 := plot(fun1, x=0..1, color = [blue], linestyle = [solid], thickness = [3], legend = ["H"]):
p2 := plot(fun2, x=0..1, color = [red], linestyle = [dash], thickness = [3], legend = ["L"]):
plots:-display({p1, p2}, title="Multiple plots");
end proc:

Explore(myPlotFunction(v), parameter={v=0..1});

Regards,

Omkar

WHen plotting f(x) and g(x) on same plot, and putting legend at bottom (default), the legends show horizontally. i.e. f(x) then g(x) on same line.

I'd like the legend to be stacked vertically, just like when the legend on the right or left, But keep it at bottom. But want to do all this in code. Not using any UI context tools or mouse.

Here is an example

restart;

f:=x->x^3-x^2+1;
g:=x->6-2*x-x^2;
the_legend:=[typeset("f(x) = ",f(x)),typeset("g(x) = ",g(x))]:
plots[setoptions](font=[TIMES,16], labelfont=[TIMES,18]):
the_title:="Plot of f(x) and g(x)":
plot([f(x),g(x)],x=-5..5,
     'gridlines',
     'color'=['red','blue'],
     'legend'=the_legend,
     'legendstyle'=['location'='bottom'],
     'axes'='normal',
     'title'=the_title,  
     'scaling'='unconstrained');

proc (x) options operator, arrow; x^3-x^2+1 end proc

proc (x) options operator, arrow; 6-2*x-x^2 end proc

#I want the above legend to be stacked vertically like the
#following one, but keep it at bottom

plot([f(x),g(x)],x=-5..5,
     'gridlines',
     'color'=['red','blue'],
     'legend'=the_legend,
     'legendstyle'=['location'='right'],
     'axes'='normal',
     'title'=the_title,  
     'scaling'='unconstrained');

 

 

Download legend_question.mw

I was wondering which theorem the following result is based on, and what the name of the sequence used is.

restart; with(LinearAlgebra)

K := proc (i::integer, j::integer) local M; M := Matrix(5, 5); M[i, i] := 1; M[j, j] := 1; return M end proc

F := proc (r, c, g, h) options operator, arrow; Adjoint(A.K(r, c).B.K(g, h).C) end proc

d1 := 3; d2 := 5; A := Matrix(d1, d2, symbol = a); B := Matrix(d2, d2, symbol = b); C := Matrix(d2, d1, symbol = c)

d1 := 3

 

d2 := 5

 

Matrix(3, 5, {(1, 1) = a[1, 1], (1, 2) = a[1, 2], (1, 3) = a[1, 3], (1, 4) = a[1, 4], (1, 5) = a[1, 5], (2, 1) = a[2, 1], (2, 2) = a[2, 2], (2, 3) = a[2, 3], (2, 4) = a[2, 4], (2, 5) = a[2, 5], (3, 1) = a[3, 1], (3, 2) = a[3, 2], (3, 3) = a[3, 3], (3, 4) = a[3, 4], (3, 5) = a[3, 5]})

 

Matrix(5, 5, {(1, 1) = b[1, 1], (1, 2) = b[1, 2], (1, 3) = b[1, 3], (1, 4) = b[1, 4], (1, 5) = b[1, 5], (2, 1) = b[2, 1], (2, 2) = b[2, 2], (2, 3) = b[2, 3], (2, 4) = b[2, 4], (2, 5) = b[2, 5], (3, 1) = b[3, 1], (3, 2) = b[3, 2], (3, 3) = b[3, 3], (3, 4) = b[3, 4], (3, 5) = b[3, 5], (4, 1) = b[4, 1], (4, 2) = b[4, 2], (4, 3) = b[4, 3], (4, 4) = b[4, 4], (4, 5) = b[4, 5], (5, 1) = b[5, 1], (5, 2) = b[5, 2], (5, 3) = b[5, 3], (5, 4) = b[5, 4], (5, 5) = b[5, 5]})

 

Matrix(%id = 36893489807004764868)

(1)

simplify(Adjoint(A.B.C)-add([F(1, 2, 1, 2), F(1, 2, 1, 3), F(1, 2, 1, 4), F(1, 2, 1, 5), F(1, 2, 2, 3), F(1, 2, 2, 4), F(1, 2, 2, 5), F(1, 2, 3, 4), F(1, 2, 3, 5), F(1, 2, 4, 5), F(1, 3, 1, 2), F(1, 3, 1, 3), F(1, 3, 1, 4), F(1, 3, 1, 5), F(1, 3, 2, 3), F(1, 3, 2, 4), F(1, 3, 2, 5), F(1, 3, 3, 4), F(1, 3, 3, 5), F(1, 3, 4, 5), F(1, 4, 1, 2), F(1, 4, 1, 3), F(1, 4, 1, 4), F(1, 4, 1, 5), F(1, 4, 2, 3), F(1, 4, 2, 4), F(1, 4, 2, 5), F(1, 4, 3, 4), F(1, 4, 3, 5), F(1, 4, 4, 5), F(1, 5, 1, 2), F(1, 5, 1, 3), F(1, 5, 1, 4), F(1, 5, 1, 5), F(1, 5, 2, 3), F(1, 5, 2, 4), F(1, 5, 2, 5), F(1, 5, 3, 4), F(1, 5, 3, 5), F(1, 5, 4, 5), F(2, 3, 1, 2), F(2, 3, 1, 3), F(2, 3, 1, 4), F(2, 3, 1, 5), F(2, 3, 2, 3), F(2, 3, 2, 4), F(2, 3, 2, 5), F(2, 3, 3, 4), F(2, 3, 3, 5), F(2, 3, 4, 5), F(2, 4, 1, 2), F(2, 4, 1, 3), F(2, 4, 1, 4), F(2, 4, 1, 5), F(2, 4, 2, 3), F(2, 4, 2, 4), F(2, 4, 2, 5), F(2, 4, 3, 4), F(2, 4, 3, 5), F(2, 4, 4, 5), F(2, 5, 1, 2), F(2, 5, 1, 3), F(2, 5, 1, 4), F(2, 5, 1, 5), F(2, 5, 2, 3), F(2, 5, 2, 4), F(2, 5, 2, 5), F(2, 5, 3, 4), F(2, 5, 3, 5), F(2, 5, 4, 5), F(3, 4, 1, 2), F(3, 4, 1, 3), F(3, 4, 1, 4), F(3, 4, 1, 5), F(3, 4, 2, 3), F(3, 4, 2, 4), F(3, 4, 2, 5), F(3, 4, 3, 4), F(3, 4, 3, 5), F(3, 4, 4, 5), F(3, 5, 1, 2), F(3, 5, 1, 3), F(3, 5, 1, 4), F(3, 5, 1, 5), F(3, 5, 2, 3), F(3, 5, 2, 4), F(3, 5, 2, 5), F(3, 5, 3, 4), F(3, 5, 3, 5), F(3, 5, 4, 5), F(4, 5, 1, 2), F(4, 5, 1, 3), F(4, 5, 1, 4), F(4, 5, 1, 5), F(4, 5, 2, 3), F(4, 5, 2, 4), F(4, 5, 2, 5), F(4, 5, 3, 4), F(4, 5, 3, 5), F(4, 5, 4, 5)]))

Matrix(%id = 36893489807027993164)

(2)

nops([F(1, 2, 1, 2), F(1, 2, 1, 3), F(1, 2, 1, 4), F(1, 2, 1, 5), F(1, 2, 2, 3), F(1, 2, 2, 4), F(1, 2, 2, 5), F(1, 2, 3, 4), F(1, 2, 3, 5), F(1, 2, 4, 5), F(1, 3, 1, 2), F(1, 3, 1, 3), F(1, 3, 1, 4), F(1, 3, 1, 5), F(1, 3, 2, 3), F(1, 3, 2, 4), F(1, 3, 2, 5), F(1, 3, 3, 4), F(1, 3, 3, 5), F(1, 3, 4, 5), F(1, 4, 1, 2), F(1, 4, 1, 3), F(1, 4, 1, 4), F(1, 4, 1, 5), F(1, 4, 2, 3), F(1, 4, 2, 4), F(1, 4, 2, 5), F(1, 4, 3, 4), F(1, 4, 3, 5), F(1, 4, 4, 5), F(1, 5, 1, 2), F(1, 5, 1, 3), F(1, 5, 1, 4), F(1, 5, 1, 5), F(1, 5, 2, 3), F(1, 5, 2, 4), F(1, 5, 2, 5), F(1, 5, 3, 4), F(1, 5, 3, 5), F(1, 5, 4, 5), F(2, 3, 1, 2), F(2, 3, 1, 3), F(2, 3, 1, 4), F(2, 3, 1, 5), F(2, 3, 2, 3), F(2, 3, 2, 4), F(2, 3, 2, 5), F(2, 3, 3, 4), F(2, 3, 3, 5), F(2, 3, 4, 5), F(2, 4, 1, 2), F(2, 4, 1, 3), F(2, 4, 1, 4), F(2, 4, 1, 5), F(2, 4, 2, 3), F(2, 4, 2, 4), F(2, 4, 2, 5), F(2, 4, 3, 4), F(2, 4, 3, 5), F(2, 4, 4, 5), F(2, 5, 1, 2), F(2, 5, 1, 3), F(2, 5, 1, 4), F(2, 5, 1, 5), F(2, 5, 2, 3), F(2, 5, 2, 4), F(2, 5, 2, 5), F(2, 5, 3, 4), F(2, 5, 3, 5), F(2, 5, 4, 5), F(3, 4, 1, 2), F(3, 4, 1, 3), F(3, 4, 1, 4), F(3, 4, 1, 5), F(3, 4, 2, 3), F(3, 4, 2, 4), F(3, 4, 2, 5), F(3, 4, 3, 4), F(3, 4, 3, 5), F(3, 4, 4, 5), F(3, 5, 1, 2), F(3, 5, 1, 3), F(3, 5, 1, 4), F(3, 5, 1, 5), F(3, 5, 2, 3), F(3, 5, 2, 4), F(3, 5, 2, 5), F(3, 5, 3, 4), F(3, 5, 3, 5), F(3, 5, 4, 5), F(4, 5, 1, 2), F(4, 5, 1, 3), F(4, 5, 1, 4), F(4, 5, 1, 5), F(4, 5, 2, 3), F(4, 5, 2, 4), F(4, 5, 2, 5), F(4, 5, 3, 4), F(4, 5, 3, 5), F(4, 5, 4, 5)])

100

(3)
 

NULL

Download Adjoint353.mw

Is this account worth keeping? I only made it to try to download Maple. I am having technical problems downloading the app using a code, and was hoping that having a MapleSoft account would stop the firewall from blocking me somehow?... not my best idea, I will admit.

So is this account worth keeping? I can't delete it for some reason anyway

Just want to know if it will help me in any way with using Maple

Given equation such as 1/A=x/A, at school we are allowed to write this as 1=x by canceling A on both side.

But in Maple, even if I tell it that A is not zero, it still refuses to simplify it and cancel A from both sides.

What could be the reson for this?

eq:= 1/A = x/A;
simplify(eq) assuming A<>0

Using Mathematica it does it:

I am not looking for workaround, I know how to force this if needed, one way could be

numer(normal((lhs-rhs)(eq)))=0

gives 1 - x = 0

My question is why Maple's simplify does not simplify it automatically? Even using simplify with size option did not. Is it just weakness in simplify or is there a subtle mathematical reason behind it which I am missing?

Maple 2025.2

Hi, I want to write a function to parse an expression into coefficients and exponents as follows:


I found a similar function (factors), but it fails in the second case.

this transforamtion including two function which i try to do, but my result is so different and even is not near i did like the author mention but i don't know how reach that outcome, the importan part is the equation 2.7

 

s1.mw

Given first order nonlinear ode which is hard to solve using standard methods, the smart dsolve sometimes uses a method where it linearizes the first order ode to a linear second order ode and solves that.

Then with that solution to the linear second order ode, it is able to find the solution of the first order ode (need to resolve the constants of integration to merge them into one, but this part is easy to do).

My question is, how and what method it uses to "linearizes by differentation"? I could not find this in any textbook I have, and not able to see how it does.

Here is an example where it uses this method to solving this first order ode

restart;

Typesetting:-Unsuppress('all'); #always do this.
Typesetting:-Settings(prime=x,'typesetprime'=true); #this says to use y'(x) instead of dy/dx    
Typesetting:-Suppress(y(x)); # this says to use y' and not y'(x)

ode:=diff(y(x),x) = (-y(x)^2+4*a*x)^2/y(x); 
infolevel[dsolve]:=5;
dsolve(ode,y(x), singsol=all);

Tracing says 

So it says, if I understand, that it differentiated the original given first order ode and then linearized the resulting second order ode to the above, which is    y''=-64 a^2 x^2 y - 16 a x y' which is certainly linear second order ode and now can be solved using kovacic algorithm. Now the solution to the first order ode can be obtained.

But when differentiating the first order ode, this is the result

expand(diff(ode,x))

So the question is, did Maple mean it "linearized" the above to y''=-64 a^2 x^2 y - 16 a x y' 

If so, then how? Did it use Taylor series? but if so, where it expanded about? Or did it use some other method?  One thing I noticed is that by multiplying the RHS above by y^2 it becomes

And "removing" the nonlinear terms (say expansion is around y=0, so higher order y terms are very small and can be removed) the RHS above becomes

    y'' y^2 = -16 y' a^2 x^2 + 32 a^2 x y

Which is close to Maple shows, but 32 instead of 64, and what about the y^2 on the left side?  Is this method even valid mathematically to do? Since solution that result will be correct only near y=0?

So far, trying to step in the debugger to find how, I was not able to find where it does that but will keep trying.

Any idea what Maple means by linearization to 2nd order and how it does it?

ps. only case I know about, where nonlinear first order ode can be transformed to linear second order ode is the Riccati ode using transformation y=u'/u. But this  first order ode is not Riccati.

3 4 5 6 7 8 9 Last Page 5 of 2451