There is another way to do this. You can create the sum of a collection of terms by converting a list (or set) to a sum as follows:

test2 := j->seq(x[j-i mod 24], i = 0 .. 8);
j -> seq(x[`mod`(j - i, 24)], i = 0 .. 8)
test2(3);
          x[3], x[2], x[1], x[0], x[23], x[22], x[21], x[20], x[19]
convert( [ test2(3) ], `+` );
      x[3] + x[2] + x[1] + x[0] + x[23] + x[22] + x[21] + x[20] + x[19]

This is an old trick that I have not had the need to use in quite some time. But, it was fun to see it return from remote storage when I read this post.

I hope this is helpful,

Doug

---------------------------------------------------------------------
Douglas B. Meade  <><
Math, USC, Columbia, SC 29208  E-mail: mailto:meade@math.sc.edu
Phone:  (803) 777-6183         URL:    http://www.math.sc.edu

Rather than just trying to understand what your MATLAB code is supposed to do and then determining how to best implement that algorithm in Maple, could you just tell us what problem you are trying to solve?

Doug

---------------------------------------------------------------------
Douglas B. Meade  <><
Math, USC, Columbia, SC 29208  E-mail: mailto:meade@math.sc.edu
Phone:  (803) 777-6183         URL:    http://www.math.sc.edu

Sure. You have h on the RHS of the assignment to h. For starters, I would change the name to which you assign the double sum. Then, I think everything goes smoothly:

H := Sum(Sum(n!/m!/j!/(n-m-j)!*(1+p/20-c-h)^m*(h-p/20)^j*(c^(n-m-j))*j*(w+r*p),j=0..n-m),m=0..n);
    n   /n - m 
  ----- |----- 
   \    | \    
    )   |  )   
   /    | /    
  ----- |----- 
  m = 0 \j = 0 

                                   m           j                         \
                 /    1           \  /    1   \   (n - m - j)            |
    factorial(n) |1 + -- p - c - h|  |h - -- p|  c            j (w + r p)|
                 \    20          /  \    20  /                          |
    ---------------------------------------------------------------------|
               factorial(m) factorial(j) factorial(n - m - j)            |
                                                                         /
value( H );
                                                             (n - 1) 
        1                          (n - 1) /-p + 20 c + 20 h\        
        -- n (20 h - p) (w + r p) c        |----------------|        
        20                                 \       c        /        

                            (n - 1)
          /       1        \       
          |----------------|       
          \-p + 20 c + 20 h/       
simplify( %, symbolic );
                          1                        
                          -- n (20 h - p) (w + r p)
                          20                       

Is this what you are hoping to get?

Doug

---------------------------------------------------------------------
Douglas B. Meade  <><
Math, USC, Columbia, SC 29208  E-mail: mailto:meade@math.sc.edu
Phone:  (803) 777-6183         URL:    http://www.math.sc.edu

I think I can see what you are trying to do. The main problem you are having is that you need to get rid of the big-O terms in the series. A simple way to do this is to convert the expression to a polynomial:

restart;
f(a) :=0: 
A := taylor(f(x), x = a, 4):
C := e[n]-A/D(f)(a):          # evaluate derivative at x=a
subs(x-a = e[n], C):
convert( %, polynom );
                          1                    2   1                    3
           D(f)(a) e[n] + - @@(D, 2)(f)(a) e[n]  + - @@(D, 3)(f)(a) e[n] 
                          2                        6                     
    e[n] - --------------------------------------------------------------
                                      D(f)(a)                            
collect( %, e[n] );
                                     3                      2
                  @@(D, 3)(f)(a) e[n]    @@(D, 2)(f)(a) e[n] 
                - -------------------- - --------------------
                       6 D(f)(a)              2 D(f)(a)      

I'm not quite sure why you are looking at C, or assuming f(a)=0. Can't you get to the same conclusion by doing?

restart;
A := taylor(f(x), x = a, 4):
C := (f(a)+D(f)(a)*e[n])-A:
subs(x-a = e[n], C);
                           /                      1                    2
     f(a) + D(f)(a) e[n] - |f(a) + D(f)(a) e[n] + - @@(D, 2)(f)(a) e[n] 
                           \                      2                     

          1                    3    /    4\\
        + - @@(D, 3)(f)(a) e[n]  + O\e[n] /|
          6                                /
convert( %, polynom );
                1                    2   1                    3
              - - @@(D, 2)(f)(a) e[n]  - - @@(D, 3)(f)(a) e[n] 
                2                        6                     

PostScript (added later): As I thought more about this as I was doing other things,
this sequence of steps is rather empty. It never makes any use of the Newton iteration.
Maybe this is what you were trying to do with the B. If so, you should be able to use
the ideas I've demonstrated to get your desired result.

Doug

---------------------------------------------------------------------
Douglas B. Meade  <><
Math, USC, Columbia, SC 29208  E-mail: mailto:meade@math.sc.edu
Phone:  (803) 777-6183         URL:    http://www.math.sc.edu

I'm not sure how your post relates to the one you quoted, but I do see at least one problem with what you have done.

You need to load the plots package once, and not in an assignment.

Try this:

with(plots);
plot1:=inequal({f,g,h,j}, phi_a=0..1, phi_Ia=0..1, optionsfeasible=(color=blue));
plot2:=iinequal({k,l,m,n}, phi_a=0..1, phi_Ia=0..1, optionsfeasible=(color=green));
plot3:=iinequal({s,t,u,v}, phi_a=0..1, phi_Ia=0..1, optionsfeasible=(color=red));
plot4:=inequal({w,x,y,z}, phi_a=0..1, phi_Ia=0..1, optionsfeasible=(color=pink));
plot5:=inequal({b,p,d,r}, phi_a=0..1, phi_Ia=0..1, optionsfeasible=(color=pink));
display([plot1, plot2, plot3, plot4, plot5]);

You won't see the individual plots as they are created, only the final one. If you want to see one of the plots, enter, e.g.,

plot1;

I hope this is helpful,

Doug

---------------------------------------------------------------------
Douglas B. Meade  <><
Math, USC, Columbia, SC 29208  E-mail: mailto:meade@math.sc.edu
Phone:  (803) 777-6183         URL:    http://www.math.sc.edu

Welcome to the perils of implicit multiplication, and "constant functions".

This has been discussed many times in MaplePrimes (and will be discussed many more times, I predict). See www.mapleprimes.com/forum/solvingvariable#comment-24852 .

To address your specific question, we have to begin by noting that you have a typo in your message Where you wrote (3+7) you most likely meant (3+4), aka 7.

Let's look at how Maple sees the two inputs:

(3+4)(7+8)

Here Maple interprets the consective parentheses (with no space) as a value of a function:  The (3+4) is a function and the (7+8) is the argument of that function. Constants are all treated as constant functions, so what Maple really sees is 7(15) which is the constant function (7) evaluated when its argument is 15. This gives the value: 7.

(3+4)*(7+8)

With the * between the parentheses Maple interprets this as multiplication: 7*15 which is 105.

Where this really gets murky is if you omit the *, and write (3+4) (7+8) (with a space between the parentheses). If your input is in 1D Maple Notation, then this is function evaluation (returning 7).  But in 2D Input the space is treated as implicit multiplication and Maple will respond with 105.

All of this is well-documented, but mistakes in these areas can be very difficult to catch. I know from experience!

I hope this is helpful,

Doug

---------------------------------------------------------------------
Douglas B. Meade  <><
Math, USC, Columbia, SC 29208  E-mail: mailto:meade@math.sc.edu
Phone:  (803) 777-6183         URL:    http://www.math.sc.edu

Nice solution.

When looking at the implementation, however, I think the construction is a little easier if the initial and terminal indices are included as "breakpoints":

breaker2:= proc(L, A, B)
        local R, n, breakpoints;
        R:= NULL;
        n:= nops(L);
        breakpoints:= [ 0,
                        select(i -> (L[i]=A and L[i+1]=B), [$1..n-1])[],
                        n ];
        if n = 2 then L
        else seq(L[breakpoints[i-1]+1..breakpoints[i]], i=2..nops(breakpoints))
        fi
    end proc;

This has not functional advantages, or disadvantages over Robert's soluiton. (Removing m is not significant, and might even make the code more difficult for some to understand.)

Doug

---------------------------------------------------------------------
Douglas B. Meade  <><
Math, USC, Columbia, SC 29208  E-mail: mailto:meade@math.sc.edu
Phone:  (803) 777-6183         URL:    http://www.math.sc.edu

When you think about what you need to do, step-by-step, this will work out very nicely.

First, your initial IVP has three parameters (F, V, and nD0). You need to tell Maple about this when you ask for a solution. Also, you want to get the "solution" as a list of procedure (you'll see why soon):

dsysT := { diff(nT(t),t)=-(F/V)*nT(t), nT(0)=nD0 };
                    /              d            F nT(t)\ 
                   { nT(0) = nD0, --- nT(t) = - ------- }
                    \              dt              V   / 
dsnT := dsolve( dsysT, numeric, parameters=[F,V,nD0], output=listprocedure );
            [t = proc(t)  ...  end;, nT(t) = proc(t)  ...  end;]

What you really want from this is the procedure that returns nT for a given numeric value of t. This is why we needed a list of procedures:

nTs := eval( nT(t), dsnT );  # you could use rhs( dsnT[2] ) but this is less clear
proc(t)  ...  end;

Now, give values to the parameters and Maple will give you a numerical approximation to the solution for a given value of t:

nTs(parameters=[1,-2,1]);
                         [F = 1., V = -2., nD0 = 1.]
nTs(2);
                             2.71828133411966455

Now, in your next IVP you can refer to the solution of the first IVP as nTs(t):

dsysD := { diff(nD(t),t)=(F/V)*(nDi-nD(t))-nD(t)*nTs(t)*sigma, nD(0)=nD0 };
       /              d          F (nDi - nD(t))                     \ 
      { nD(0) = nD0, --- nD(t) = --------------- - nD(t) nTs(t) sigma }
       \              dt                V                            / 
dsnD := dsolve(dsysD, nD(t), numeric, parameters=[F,V,nDi,sigma,nD0], output=listprocedure );
            [t = proc(t)  ...  end;, nD(t) = proc(t)  ...  end;]
nDs := eval( nD(t), dsnD );
proc(t)  ...  end;
nDs( parameters=[1,-2,3,4,5] );
              [F = 1., V = -2., nDi = 3., sigma = 4., nD0 = 5.]
nDs( 1 );
                            -0.223125933665962451

Note that the parameter values for the T IVP must be specified separately from those for the D IVP (they could even be different - Maple won't check for this).

Now you should be able to generalize this for your real problem. (If you want to refer to a derivative of a solution, either eliminate the derivative by using the RHS of the appropriate ODE or extend the problem to a system in which the derivative is a new dependent variable and find an ODE for this variable.)

I hope this helps,

Doug

---------------------------------------------------------------------
Douglas B. Meade  <><
Math, USC, Columbia, SC 29208  E-mail: mailto:meade@math.sc.edu
Phone:  (803) 777-6183         URL:    http://www.math.sc.edu

If you want to plot a collection of explicit points, you might want to look at pointplot (in the plots package).

with( plots ):
pointplot( [ [1,2], [3,4] ], axes=boxed );

Or, just do it with regular plot:

plot( [ [1,2],[3,4] ], style=point, axes=boxed );

See the online help for either for more information, and examples.

Doug

---------------------------------------------------------------------
Douglas B. Meade  <><
Math, USC, Columbia, SC 29208  E-mail: mailto:meade@math.sc.edu
Phone:  (803) 777-6183         URL:    http://www.math.sc.edu

It is good advise to use Matrix rather than matrix.

The command you can use to "polish" your matrix is fnormal.

with( LinearAlgebra ):
A:=Matrix(2,2,[4, 12, 2, 0.04]);
map(fnormal,A,3);

The way fnormal( e, d ) works is that it treats all numbers in e that are less than 10^(-d+2) as if they are zero. In your case,  your threshold was 0.1=10^(-1), so I use d=3. For cutoffs that are not powers of 10, you might want to specify the third optional argument to fnormal. (See ?fnormal for more details - but no example showing how this works.)

There is also a Map command that can be useful when working with Matrices. It's very similar to map, except that it works "in place". That is, the results replace the original contents of the Matrix. Try, for example,

A;
Map( fnormal, A, 3 );
A;

Notice that the Matrix A has been changed by the Map command.

Doug

---------------------------------------------------------------------
Douglas B. Meade  <><
Math, USC, Columbia, SC 29208  E-mail: mailto:meade@math.sc.edu
Phone:  (803) 777-6183         URL:    http://www.math.sc.edu

This effect is not difficult to reproduce in Maple (up to a point).

First, load the plots package.

with( plots ):

Next, assign a name to each plot you create:

p1 := plot( sin(x), x=0..Pi );

To see the plot:

p1;

Create the next part of the graph you are creating:

p2 := plot( sin(2*x), x=0..Pi, color=blue );

To create a composite plot:

display( p1, p2 );

And so on,

p3 := plot( sin(3*x), x=0..Pi, color=green );
display( p1, p2, p3 );

I'll close by noting that you can also create an animation in the same manner:

display( [p1,p2,p3], insequence=true );

This is not quite the same as MATLAB's hold on, but it give much of the same functionality (and sometimes more!).

Doug

---------------------------------------------------------------------
Douglas B. Meade  <><
Math, USC, Columbia, SC 29208  E-mail: mailto:meade@math.sc.edu
Phone:  (803) 777-6183         URL:    http://www.math.sc.edu

Here is how I would continue to do this, with animate. Define your own command to create individual frames in the animation:

PlotFrame := i -> plots:-pointplot3d( plotter[i], style=line, linestyle=solid, thickness=5 );

Once you have debugged PlotFrame, you can call animate as follows:

animate( PlotFrame, [c], c=[$1..N], trace=N-1 );

where N is the total number of frames. Note the use of a list to avoid the problem Robert pointed out with Maple using a float for the parameter range. (Maybe you can replace N with an appropriate dimension of plotter.)

Lastly, the trace option should be set to N-1 (not 1 as I wrote in my original post). The total number of frames to be shown at the end of the animation is one more than this setting, so N-1 means that all N frames will be shown at the end.

I think this should be pretty close to what you are trying to do.

Doug

---------------------------------------------------------------------
Douglas B. Meade  <><
Math, USC, Columbia, SC 29208  E-mail: mailto:meade@math.sc.edu
Phone:  (803) 777-6183         URL:    http://www.math.sc.edu

The op command returns the argument(s) of a function. In your case, you want

op( sol );

If there are multiple arguments, op returns all arguments as an expression sequence. To get the 2nd argument of an object, use op(2, ... ).

Doug

---------------------------------------------------------------------
Douglas B. Meade  <><
Math, USC, Columbia, SC 29208  E-mail: mailto:meade@math.sc.edu
Phone:  (803) 777-6183         URL:    http://www.math.sc.edu

 animate3d has been replaced with animate. The animate command is pretty powerful but can take a little time to learn to use.

Suppose you generate the individual frames with the command

pointplot3d( plotter[1], other options );

then you would create the animation of the 20 frames with

animate( pointplot3d, [plotter[i], other options], i=1..20 );

Now, it sounds as though you want each frame to simply add one more set of points to the overall image. For this, use the trace option of the animate command:

animate( pointplot3d, [plotter[i], other options], i=1..20, trace=1 );

There is also a background option that allows you to specify a common background for all frames in the animation.

See the help for plots,animate for complete details and some examples.

Doug

---------------------------------------------------------------------
Douglas B. Meade  <><
Math, USC, Columbia, SC 29208  E-mail: mailto:meade@math.sc.edu
Phone:  (803) 777-6183         URL:    http://www.math.sc.edu

Can you make available this worksheet? That would allow us to try different things and make specific suggestions.

How are  your arrays created? It might not be very difficult to convert these to an Array or Matrix.

Based on what I am seeing here, isn't xcoord[1] simply the number 3? What if you try the following?

with(plots):
pointplot3d( xcoord, ycoord, zcoord );

Doug

---------------------------------------------------------------------
Douglas B. Meade  <><
Math, USC, Columbia, SC 29208  E-mail: mailto:meade@math.sc.edu
Phone:  (803) 777-6183         URL:    http://www.math.sc.edu