We can use sscanf to break the input string into word-sized chunks, and then pass those words to sprintf so that they can be reassembled with the spaces between words. The syntax of sprintf and printf is identical. The reason that I chose sprintf is that there's no way to programmatically access the output of printf; it produces output for viewing only. 

restart;
S:= "IDIDNTSEETHATDOUBLERAINBOW":
wordlength:= 5:
# Construct scan format code
WFI:= sprintf("%%%ds", wordlength);

                         "%5s"

# Construct print format code
WFO:= sprintf("%%0.%ds", wordlength);

                          "%0.5s"

# Count words
words:= iquo(length(S), wordlength, 'r'):
if r>0 then  words:= words+1  end if:

sprintf(cat(WFO, cat(" ",WFO) $ words-1), sscanf(S, cat(WFI $ words))[]);
                            
               "IDIDN TSEET HATDO UBLER AINBO W"

Is this what you were looking for?

The error message about "points with fewer or more than 2 components" clearly provides a clue to the problem. The plot command is two-dimensional (2D), and thus cannot handle the 3D points that you try to give it in the plots axyz and ixyz. For these two plots, change the command plot to pointplot3d.

There's another bug in your code, more subtle and insidious. You use variable n several ways, some of which are in conflict. You use it as a parameter in the first differential equation de1, and as such it also appears in the f:= (x,y,z)-> ... in procedures ESys and Fourth. But you also use it as a local variable (for a loop index) in those procedures. You need to change n to something else, say nn, in four places: the definition of de1, the line where you set the numeric values of the parameters, and the two f:= (x,y,z)-> ... lines.

Consider adding some spaces to your code to make it easier to read.

If you have further questions about this piece of code, please upload a worksheet. This piece is a bit too long for an easy cut-and-paste by the reader (me in this case). After I cut-and-paste, I still have to manually block off the comments and remove the extra >s.

10^5 is not too many points for the plot renderer, but 10^6 is unwieldy. Execute these commands to create the plot:

nn:= nextprime(10^7):
zz:= 1:
N:= 10^5:
Z:= Matrix(N, 2):
for k to N do
   zz:= (zz^2+1) mod nn;
   Z[k,1]:= k:  Z[k,2]:= zz
end do:
plot(Z, style= point);

The plot will appear as a solid block of points. Put the mouse pointer in the plot and right click. Select Manipulator from the context menu, then Scale from the submenu. Select a rectangle on the plot which is horizontally narrow but vertically covers the full extent of the plot. Do this a few more times, until you're satisfied with the narrowness of the horizontal range. Then go to the manipuator submenu again and select Pan. Now the plot is horizontally scrollable.

The problem is that EmpiricalDistribution is a discrete distribution, yet the sample is drawn from a continuous distribution. Something almost analogous to EmpiricalDistribution that works for continuous distributions is KernelDensity. Specifically, KernelDensity returns a probability density function (PDF) deduced from a sample drawn from a continuous distribution. Once you have the PDF, you can create a custom Distribution, and then a RandomVariable from that.

Download KernelDensity.mw

Answering your first question, I strongly suspect that in the code before the code that you show that you assigned a plot to the variable a. Doing a restart should fix that.

I believe that the following will work in Maple 12. Let me know if it does not.

Download Euler.mw

Here (attached) is your code with my (minor) modifications to make it run in hardware floating point. I ran it up to size 512 in less than a minute, and got reasonable accuracy (norm of the residuals 1e-16 seems reasonable to me). Let me know if this fulfills your needs.

Markov_a2.mw

In your program, what is the purpose of the following lines?

if type(y, rational) = false then
   y:= evalf(y);
   y:= convert(y, rational)
end if;

In other words, why are you using exact rational arithmetic, even after floating point approximation?

I ran your code for m=50, and the answers have more than 10000 digits in the numerator and denominator. It could be over a million digits for m=300.

Maple's sign is the sign of the leading coefficient of a polynomial, so sign(x) is always 1. Maple's signum is fine if you are dealing with strictly real expressions; for complex values, it may be counterintuitive (I'm not saying that it's wrong---it just may not be what you expect). An alternative is csgn. See ?signum, ?csgn, and ?sign.

Your IVP can be solved symbolically. The solution is x(t) = (1+exp(-t))/2. Thus .5 < x(t) < 1 for all t > 0. Thus both piecewises and the abs make no difference at all and the whole ODE reduces to

diff(x(t), t) = .5 - x(t).

When you make assumptions on a variable, a copy of that variable is made. If you then assign to that same variable, then you'll have the original and the copy being used in different ways, and things get confusing. The takeaway advice for the inexperienced user is Don't assign to assumed variables.

Your last two lines were

• a:=solve(z1,a);
• z2;

You should change those to

• soln:= solve(z1, {a});
• eval(z2, soln);

When used in this way, eval has most of the desired effect of an assignment statement without the side-effect of changing the global state. You can of course assign the result of eval to another variable.

Better than what you asked for---the solutions for a static list---we can make a "blackbox" solution that provides the x for any given y. This is what symbolic computation systems excel at.

Download blackbox.mw

A solution is to ?overload the procedure `diff/Int`.

I coded the overload of `diff/Int` solution. It was actually much easier than I was anticipating, just 11 lines of code. This solution illustrates a very generally applicable technique, one that can be used to modify the behavior of almost any Maple library or even built-in procedure. This technique is not often discussed on MaplePrimes, so I suggest that most readers save this Answer in their Favourites.

restart;

`diff/Int/orig`:= eval(`diff/Int`):
unprotect(`diff/Int`);
`diff/Int`:= overload([
    proc(f::algebraic, x::{name, name=range}, M::(identical(method)=name))
        option overload;
        subsindets(
            diff(Int(f,x), args[-1]),
            specfunc(anything, 'Int'),
            subs([_DUMMY= args[3..-2]], J-> 'Int'(op(J), _DUMMY))
        )
    end proc,
    
    `diff/Int/orig`
]):

Test

I have only overloaded the case where the first argument to Int is type algebraic, the second argument is a name or name= range, and the third argument is method= name. All other cases should behave as they did originally. If you need the overload to work for more general cases of Int, please let me know.

Download diffInt_overload.mw

You wrote:

Notice that the order of the output from the EulerLagrange command when three functions are input does not correspond to the order of the functions listed, and when they are requested separately, the last function is not included:

The two phenomena are unrelated. As to the order: The EulerLagrange command always gives its result as a set, in curly braces. There is no way to control the order of elements in a Maple set: An internal lexicographic ordering is used so that set-membership checks can be done with binary search and set order is preserved across sessions. (Older Maple used an address-based ordering and set order varied by session.)

You wrote:

T3:=(1/2)*m*(diff(r(t),t)^2+r(t)^2*diff(theta(t),t)^2+diff(z(t),t)^2):
U3:=m*g*z(t):
L3:=T3-U3;

As for the second issue: Your equations above use diff(theta(t), t) but they never use theta(t) undifferentiated. This confuses EulerLagrange, and perhaps it should be rejected as invalid input. If you change diff(theta(t),t) to theta_prime(t) and do not include this in the solution variables, then the union of the individual results for z(t) and r(t) will be the same as the result when both of those are requested.

It can be done like this (I just picked a simple curve for the implicitplot).

P:= plots:-implicitplot(y= exp(x), x= -10..10, y= 0..10^4):
data:= op([1,1], P):
plots:-logplot(convert(data, listlist));