uses the code

  restart;
  ineqs1:= [ y > 0,
             y < 1/2*(2*x+x^2)+1/2*sqrt(4*x^3+x^4),
             x > 0,
             x < 1/13*(5-2*sqrt(3))
           ]:
  p1:= plots:-inequal( ineqs1,
                       x=0..1,
                       y=0..1,
                       nolines,
                       optionsimplicit=[grid=[100,100]],
                       color=red,
                       transparency=0.5
                     ):
  ineqs2:= [ y > 1/2*(2*x+x^2)+1/2*sqrt(4*x^3+x^4),
             y < 2*x,
             x > 0,
             x < 1/13*(5-2*sqrt(3))
           ]:
  p2:= plots:-inequal( ineqs2,
                       x=0..1,
                       y=0..1,
                       nolines,
                       optionsimplicit=[grid=[100,100]],
                       color=blue,
                       transparency=0.5
                   ):
  r1:=op~(plottools:-getdata(p1, rangesonly)):
  r2:=op~(plottools:-getdata(p2, rangesonly)):
  plots:-display([p1,p2],
                  view=[ floor( 100*min( r1[1], r2[1]))/100..ceil( 100*max( r1[2], r2[2]))/100,
                         floor( 100*min( r1[3], r2[3]))/100..ceil( 100*max( r1[4], r2[4]))/100
                       ],
                  legend=[ "ineqs1", "ineqs2"],
                  legendstyle=[font=[times, roman, 16]]
                );

to produce

See the attached worksheet

ineqplt.mw

Obviously my selection of various plotting options in the above is pretty arbitrary - feel free to change them!

As Carl has noted, your code doesn't work, because it requires 'recursive' assignments

If I make some guesses about what you actually intended, and fix some other inefficiencies (ie whenver you are adding a finiet number of terms, use add(), not sum(), etc,etc), then I can make the code in the attached run.

The output from the attached worksheet can be summarised as

N       Computes            cputime             bytesused
1             u[1]                  0.063               2,027,208
2         u[1]..u[2]              0.031               4,328,664
3         u[1]..u[3]              1.763           290,992,944
4         u[1]..u[4]          624.753      39,235,739,664

I didn't calculate any further, because my life is too short!!!!

Download usage.mw

why you did not use

restart;
f:= sqrt(x);    # assumed to be unknown
u:= convert(f, string);
s:= StringTools:-Encode( u, ':-base64');
v:= StringTools:-Decode( s, 'encoding' = ':-base64' );
parse(v);

??? Or am I missing something

The above code produces

It is in fact possible to find the equations of motion relatively easily, without recourse to Maple's built-in Euler-Lagrange() command or the Physics package - although obviously there is nothing wrong with either of these approaches. (I'd probably use the Euler-Lagrange() command for preference!).

One only has to

1. make some variable name substitutions to facilitate computing the required partial derivatives
2. reverse the variable name substitutions
3. compute the required and set up the Euler-Lagrange equations

Code for this is as shown

  restart;
#
# The original Lagrangian. NB using the code provided
# in Rouben Rostamian original response
#
  L := J/2*(diff(x(t),t)/R)^2
       + m__r/2*diff(x(t),t)^2
       + m__p/2*((diff(x(t),t) + l*diff(phi(t),t)*cos(phi(t)))^2 + (l*diff(phi(t),t)*sin(phi(t)))^2)
       + m__p*g*l*cos(phi(t));
#
# Define a variable substitution which facilitates the
# generation of the Euler-Lagrange equations.
#
# Note that the 'order' of these substutions is important
# one has to make sure that the derivative terms are done
# 'first'
#
  subVars:= [ diff(x(t),t)=xdot,
              diff(phi(t), t)=phidot,
              x(t)=x,
              phi(t)=phi
            ]:
# 
# Substitute these in the above Lagrangian
#
  L:= subs( subVars,
            L
          );
#
# Define the reverse substitution
#
  backSubs:= rhs~(subVars)=~(lhs~(subVars)):
#
# Compute the Euler-Lagrange equations.
# This comprises
#
# 1. Compute diff(L, x), diff(L, phi)
#    diff(L, xdot), diff(L, phidot)
# 2. Make the reverse substituions defined
#    in backSubs above
# 3. Then compute the EL equations using
#
#    diff( diff(L, xdot), t)=diff(L, x)
#    diff( diff(L, phidot), t) =diff(L, phi)
#
   EL1:= diff(subs( backSubs, diff(L, xdot)),t) = subs( backSubs, diff(L, x));
   EL2:= diff(subs( backSubs, diff(L, phidot)),t) = subs( backSubs, diff(L, phi));

The attached worksheet produces the same answers as Rouben Rousamian obtained.

Download lagrange.mw

which way round yuwant to do the composition - but the following code does it both ways, so just pick the one you want

  restart;
  with(PolynomialTools):

  g:= FromCoefficientList([1,2,3],x);
  h:= FromCoefficientList([4,5,6,7],x);

  expand( subs( x=g, h) );
  CoefficientList
  ( subs
    ( x=g, h),
    x
  );

  expand( subs( x=h, g) );
  CoefficientList
  ( subs
    ( x=h, g),
    x
  );

the following possibilities should work - depending on whether you are actually turning off 'labels' and/or 'tickmarks'

  restart;
#
# The default plot
#
  p1:=plot(sin(x), x=0..2*Pi );
#
# Turn off the axes labels
#
  p2:=plot(sin(x), x=0..2*Pi, labels=["", ""]);
#
# Turn off the axes labels and the tickmarks
#
  p3:=plot(sin(x), x=0..2*Pi, labels=["", ""], axis=[tickmarks=0]);

a sixth-order polynomial doesn't look too bad,

The code

  restart;
  with(Statistics):
  Y := Vector([1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 2, 2, 3, 8, 12, 12, 22, 30, 40, 44, 51, 65, 70, 97, 111, 131, 135, 174, 184, 210, 214, 232, 238, 254, 276, 285, 305, 318, 323, 343, 373, 407, 442, 493, 542, 627, 665, 782, 873, 981, 1095, 1182, 1273, 1337, 1532, 1728, 1932, 2170, 2388, 2558, 2802, 2950, 3145, 3526, 3912, 4151, 4399, 4641, 4787, 4971, 5162, 5445, 5621, 5959, 6175, 6401, 6677, 7016, 7261, 7526, 7839, 8068, 8344, 8733, 8915, 9302, 9855, 10162, 10819, 11166, 11516, 11844, 12233, 12486, 12801, 13464, 13873, 14554, 15181, 15682, 16085, 16658, 17148, 17735, 18480, 19147, 19808, 20244, 20919, 21371, 22420, 22614, 23298, 24077, 24567, 25133, 25694, 26484, 27110, 27564, 28162, 28711, 29268, 29789, 30249, 30784, 31323, 31987, 32558, 33153, 33616, 34259, 34854, 35454, 36107, 36663, 37225, 37801, 38344, 38948, 39539, 39977, 40532, 41180, 41804, 42208, 42689, 43151, 43537, 43841, 44129, 44433, 44890, 45244, 45687, 46140, 46577, 46867, 47290, 47743, 48116, 48445, 48770, 49068, 49485, 49895, 50488, 50964, 51304, 51905, 52227, 52584, 52800, 53021, 53317, 53477, 53727, 53865, 54008, 54247, 54463, 54588, 54743, 54905, 55005, 55160, 55456, 55632, 55829, 56017, 56177, 56256, 56388, 56478, 56604, 56735, 56956, 57145, 57243, 57437, 57613, 57724, 57849, 58062, 58198, 58324, 58460, 58647, 58848, 59001, 59127, 59287, 59345, 59465, 59583, 59738, 59841, 59992, 60103, 60266, 60430, 60655, 60834, 60982, 61194, 61307, 61440, 61558, 61630, 61667, 61805, 61882, 61930, 61992, 62111, 62224, 62371, 62521, 62691, 62853, 62964, 63036, 63173, 63328, 63508, 63731, 63790, 64090, 64184, 64336, 64516, 64728, 64884, 64996, 65148, 65305, 65693, 65693, 65839, 65982, 66228, 66230, 66439, 66607, 66805, 66974, 67220, 67330, 67412, 67557, 67838, 67960, 68303, 68627, 68937, 69225, 69645, 70195, 70609, 71344, 72140, 72757, 73175, 73394, 74132, 75062, 76207, 77013, 77933, 78434, 78790, 79789]):
  polyOrder:=6:
  f1:=Fit(add( a[j]*t^j,j=0..polyOrder), <seq(j, j=1..numelems(Y))>, Y, t, summarize=true);
  plots:-display( [ plot( <seq(j, j=1..numelems(Y))>, Y, color=blue),
                    plot( f1, t=1..numelems(Y), color=red)
                  ]
                );
  meanResid:= Statistics:-Mean( [seq( sqrt((eval(f1, t=j)-Y[j])^2), j=1..numelems(Y))]);

will produce the output

See the attached worksheet

polyFit.mw

feel free to try higher or lower values for the variable polyOrder - but my calibrated eyeball says 'order 6' is about right

Why start with the C code?

Since the Hilbert curve has a Lindemayer representation, it is relatively simple within Maple to generate, the Hilbert curve of any order.

The following code

  restart;
#
# Procedure to generate the Lindenmayer representation
# for the Hilbert curve
#
  Hcurve:= proc(n)
                uses Fractals:-LSystem, StringTools;
                local state:="A",
                      rules:=["A"=" +BF-AFA-FB+", "B"="-AF+BFB+FA-"],
                      cons:=["F"="draw:1","+"="turn:-90", "-"="turn:90"],
                      opstate;
                opstate:= SubstituteAll
                          ( SubstituteAll
                            ( Iterate(state, rules, n),
                              "A",
                              ""
                            ),
                            "B",
                            ""
                          );
                return opstate, cons
         end proc:
#
# Example - plot the Hilbert curve for 5 iterations
#
  Fractals:-LSystem:-LSystemPlot(Hcurve(5));

will  produce this plot

see the attached worksheet

Hcurve.mw

A warning - this curver really does "fill space" - and for more than about 9 iterations the returned plot is so "dense" that it appear almost ike a competely black square!!!

1. In the first case I agree with Kitonum, that the use of the solve() command to 'isolate' a specific term of an ode, is a bit 'unusual', so that I am mildly surprised that it  generally "works". If I ever wanted to do this, in what I would consider to be a 'safe' way I can see quite a lot of freeze() and thaw() commands being used!
2. However the error which user nm is reporting is not really related  to this practice. It is trivial to substitute other simple names for y(x), D(y)(x) etc, apply solve, and generate the same error.

For example the following code

restart;
interface(version);
Physics:-Version();
interface(warnlevel=4):
kernelopts('assertlevel'=2):
eq1:=exp(x)*sin(y(x))-3*x^2+(exp(x)*cos(y(x))+(1/3)/y(x)^(2/3))*(D(y))(x) = 0;
eq2:=subs( [y(x)=var1, D(y)(x)=var2], eq1);
solve(eq2, var1);

attempts to solve 'eq2', for the simple name 'var1' and generates the same error. The output is

So it would seem that there is some kind of 'bug' somewhere. I could accept that the command 'solve(eq2, var1)' might *fail* to come up with a solution, but it shouldn't produce an error!

Kitonumm's suggestions are not "wrong", but

1. the first does not typeset the final fraction correctly
2. the second essentially converts the character '3' in the final denominator to a name, which therefore prints in an italic font

The following code avoids both of these issues

restart;
P := plot(sin(x)-x+(x^3)/3!, x=-2*Pi..2*Pi) :
L := plots[textplot]( [ 3 , 17 , sin(x)-x+(x^3)/''3!''],  axes = none) :
plots[display](P ,L, axes = normal);

The Maple toolbar contains the following two (adjacent) icons for indenting and outdenting code. Use them!

you can use the 'iterations, option. Essentially this will perform the desired operation, multiple times (as defined by 'iterations'), then divide the results by the number of iterations. (I think(?) that this means you may get the "wrong" answer from the memory used!)

Compare the results from

SOL:=[1,2,3,4];
CodeTools:-Usage(SOL[1]);
CodeTools:-Usage(SOL[1], iterations=10^6);

which give

memory used=1.06KiB, alloc change=0 bytes, cpu time=0ns, real time=0ns, gc time=0ns

and

memory used=32 bytes, alloc change=0 bytes, cpu time=5.5e+02ns, real time=5.1e+02ns, gc time=78ns

respectively.

on how your code has been written, but the CodeTools:-Usage() command ought to be helpful

Difficult to work out exactly what you are trying to achieve, but one problem  you have is in the definition of anonymous procedures. I think you ought to be using unapply(expression, argument), so the expression will be evaluated before the procedure is created.

See the attached - where in the procedure Test1() I have changed the way all of the anonymous procedures are defined. This code now "runs", although whether it provides the answer you want, I have no idea

The procedure Test1(), contains a number of redundancies including a whole 'if' condition. I have removed (most?) of these redundancies in the procedure Test2() - which generates the same answer as Test1, using less code

tproc.mw

with your code.

Would someone else (like me) do it slightly differently? Yes - maybe. I'd (probably?) do ir as in the attached, where the main difference is that I use multiple 'return' statements. This (sometimes) avoids unnecessary computation associated with the use of just returning the final evaluated statement, whilst emphasing the fact that the procedure may 'exit' at various points.

Is mine 'better'? That is a subjective assessment that I wouldn't want to get into. I just find my version somewhat more 'readable'

(And due to the incompetence of the idiots who maintain this website, this simple worksheet will not display here, so you will have to download)

Download aproc.mw