Here's a procedure that does what you want:

YourProc:= proc(
     fx::algebraic,
     `x= a..b`::And(name=range(realcons), satisfies(r-> op([2,1], r) < op([2,2], r))),
     n::posint,
     {functionoptions::list({symbol, name= anything}):= []},
     {pointoptions::list({symbol, name= anything}):= []}
)
local
     x:= op(1, `x= a..b`),
     a:= op([2,1], `x= a..b`),
     r:= op([2,2], `x= a..b`) - a,   #b-a
     Df:= unapply(diff(fx, x), x),
     k,   #index
     xs:= [seq(a+k*r/n, k= 0..n)]
;
     plots:-display(
          plot(fx, `x= a..b`, functionoptions[]),
          plot([seq([k, Df(k)], k= xs)], style= point, pointoptions[]),
          _rest
     )
end proc:

And an example of its use:

YourProc(
     x^2, x= -2..2, 5,
     functionoptions= [color= green],
     pointoptions= [symbol= cross, symbolsize= 24],
     axes= boxed
);

The last three arguments are optional. You could just as well use it as

YourProc(x^2, x= -2..2, 5);

Note that intervals in Maple are specified as 2..5 rather than [2,5]. However, if you insist on using [a,b], I can easily change the procedure to accomodate that.

As a workaround, you can use Maple Input (aka 1D input).

The direct cause of your problem is that you're not allowed to assign to a parameter, in this case M. You can make a local copy of M like this:

M1:= copy(M);

Then make all assignments to M1 and return M1.

There are also a few places where you used = that should've been :=.

You should begin your procedure with 

uses LinearAlgebra;

or

uses LA= LinearAlgebra;

Don't rely on the with statement being active in your procedure---sometimes it may be, sometimes it may not be.

Suppose that your dsolve result is assigned to ans:

ans:= dsolve(...);

Then do

u:= unapply(eval(u(t), ans), t);

If g=0, then f can be arbitrary and the jacobian will still be zero. That's why the fourth solution shows f unchanged.

plot3d([r, theta, Pi/6], theta= 0..2*Pi, r= 0..1, coords= spherical);

Acer's Answer is very complete. Here's a quick-and-dirty version that works in any Maple GUI. Define two procedures to be used as binary operators:

`&<`:= (a,b)-> a*exp(I*b*Pi/180):

`&+`:= proc(a,b)
local x,y;
     (x,y):= op(evalf(polar(a+b)));
     '`&<`'(x, evalf(y*180/Pi))
end proc:

Now,

(5 &< 30) &+ (4 &< 60);

Your matrix is only Hermitian if all the variables are real. If you do substitute real values for the variables, then the imaginary parts of those two eigenvalues will be zero. Is it possible to prove, in general, that those imaginary parts are real (other than by saying that the eigenvalues came from a Hermitian matrix)? I doubt it. It's well known that there are real algebraic numbers that can't be expressed in radical form without using imaginary numbers. For example, consider the roots of x^3 - 3*x + 1.

In the worksheet below, I approximate your infinite sum using an asymptotic series for BesselJ. I just did the sum itself; I ignored the constant multipliers that you have in front of the sum.

I'm not totally confident of the mathematical validity of what I did. In particular, I have some skepticism about simplifying an asymptotic series under the assumption that the variable is integer. Hopefully, someone else can comment.

P.S.: Independent computation in later Replies confirms the accuracy of this approach, although it turns out that there's a much easier way.

Download Bessel_infinite_sum.mw

In Maple, Pi is with a capital P.

I can't conceive of any message-encoding application where you wouldn't want the same number of bits in each byte. You can get the same number with the bits option to Split:

message:= `67A`:
bitP:= Bits:-Split~(convert(message, bytes), bits= 8);

Then the total number of bits is simply

8*nops(bitP);

     24

If you insist on using a variable number of bits per byte, the total number of bits can be counted as

bitP:= Bits:-Split~(convert(message, bytes));

`+`(nops~(bitP)[]);

19

L:= [2,4,9,15,20]:
parse(cat(L[]));

@fbackelj For what it's worth, ConditionNumber(B,2) (i.e., using the 2-norm) returns the correct result:

LinearAlgebra:-ConditionNumber(B,2);

     3.225504927

Compare with

`/`(LinearAlgebra:-SingularValues(B, output= list)[[1,-1]][]);

     3.22550492667769

Regarding the condition number in the infinity-norm: Is there any way to compute the exact answer without inverting the matrix? And if one has to invert the matrix, then there's no practical reason to compute the condition number. Floating-point matrices are all about practical computation. If you want exact answers, then stick with rational matrices, or use

`*`(LinearAlgebra:-Norm~([B,1/B], infinity)[]);

      6.

and hope that the inverse isn't incorrect due to ill-conditioning.

For these reasons, I'm reluctant to call your issue a bug, although the help page should be made clearer. Perhaps an option attemptexact can be added to ConditionNumber, and it'll do the equivalent of my last line of code.

Try typing Ctrl-T to switch to text mode.

You can get a functional numeric solution like this:

SolNum:= a-> fsolve(-x^3+a*x^2-ln(x), x);

With this, you can do things like

SolNum(.1);

     0.7225317106

or

plot(SolNum, -2..2);

You can get a functional symbolic solution like this:

SolSym:= unapply(solve(-x^3+a*x^2-ln(x), x), a);

The result is not much to look at, but you can do with it symbolic things such as taking its derivative:

diff(SolSym(a), a);