To clarify Jacques' warning about positioning, consider the ODE

sys := {diff(x(t),t)=sin(t),diff(y(t),t)=cos(t),x(0)=1,y(0)=1}:

and its numeric solution

dsol := dsolve(sys, numeric):
dsol(1); 
          [t = 1., x(t) = 1.45969768100667418, y(t) = 1.84147131233092276]

You do not want to assume that x(t) will always be the second element in the list (for a particular solution, the order won't change, but a subsequent call to dsolve could conceivably change the order of x(t) and y(t), depending on how the ordering was created). To create a function of x(t) and y(t) you might do

S := tt -> eval(3*x(t) + 4*y(t), dsol(tt)):

The usual way to plot an output of dsolve/numeric is to use plots[odeplot]:

plots[odeplot](dsol, [t, 3*x(t)+4*y(t)], 0..1)

If I understand correctly, there isn't a way to do what you want. That is, given an existing module there is no straightforward way to add an export to that module. Even modifying an export of a module isn't generally possible (there are some hackish ways, but they have limits). This is a drawback of the modular approach compared to the old table-based package approach used in Maple R5—while the modular approach is more robust and better from the programmar's (creator's) viewpoint, it is less malleable for the user. It is, of course, possible to add a new module to an existing archive, though usually it is better to just create a new archive and insert it in the libname path. Use LibraryTools[Save]. To make a module ``withable,'' give it the package option. For example:

mymodule := module()
export A,B:
option package;
  A := proc() ... end proc:
  B := proc() ... end proc:
end module:
LibraryTools:-Save(mymodule,"mylib.mla"):

The ability to rebind operators, such as +, was introduced in Maple 6; however, the overload option (and procedure) was introduced in Maple 9.5.

See the help page for NextAfter.

> Digits := 20:             
> epsilon := NextAfter(0,1);
                                                                   -20
                                                  epsilon := 0.1 10

> 1+epsilon;                
                                                 1.0000000000000000000

That isn't quite what you want, but...

I'm not sure I understand, but the following (rather ugly) code does what I think you want

restart; with(stats); Digits := 14;

N := 4:
L := 0.002:
beta:=0.85:
cost := Vector([200,100,700,1000]):
lambda := Vector([0.2,0.6,0.3,0.7]):
Lambda := add(lambda[i], i=1..N):

DFR := s -> ( add(lambda[i]*add((L*lambda[i])^k*exp(-L*lambda[i])/k!
                                , k = 0 .. s[i]-1)
                  , i = 1 .. N)
              / Lambda );

IN := Matrix(N,N,shape=identity):
S := Vector[row](4):
maxval := DFR(S);

while maxval < beta do
    prevval := maxval;
    maxrelval := 0;
    for k to N do
        Sk := S + IN[k,1..N];
        val := DFR(Sk);
        relval := (val - prevval)/cost[k];
        if relval > maxrelval then
            maxrelval := relval;
            maxval := val;
            Snxt := Sk;
        end if;
    end do;
    if maxrelval = 0 then
        print("cannot increase");
        break;
    end if;
    S := Snxt;
end do:

S, maxval;

I think the following does what you want.

restart; with(stats); Digits := 14;

N := 4:
L := 0.002:
beta:=0.85:
lambda := Vector([0.2,0.6,0.3,0.7]):
Lambda := add(lambda[i], i=1..N):

DFR := s -> ( add(lambda[i]*add((L*lambda[i])^k*exp(-L*lambda[i])/k!
                                , k = 0 .. s[i]-1)
                  , i = 1 .. N)
              / Lambda );

IN := Matrix(N,N,shape=identity):
S := Vector[row](4):
maxval := DFR(S);

while maxval < beta do
    found := false; # flag to avoid infinite loop if no Sk is better than S.
    for k to N do
        Sk := S + IN[k,1..N];
        val := DFR(Sk);
        if val > maxval then
            maxval := val;
            Snxt := Sk;
            found := true;
        end if;
    end do;
    if not found then break; end if; # may want to add warning
    S := Snxt;
end do;

maxval;

I don't understand what you are attempting, however, there seem to be two clear problems. First, you probably want to use add instead of sum in the computation of DFR (you used one of each). Second, your do loop is confusing. It's index is i, but you also use i as the index in the adds. These are different i's and I'm not sure which you intend. Change one of them to a j to make it clear what you are doing.

A good method for creating types is to use the TypeTools module, in particular, its AddType export. You'd do something like:

Automatism := module()
export zeros;
local ModuleLoad,ModuleUnload;

  ModuleLoad := proc()
    TypeTools:-AddType( 'TransferFunction', handler );
  end proc;

  ModuleUnload := proc()
    TypeTools:-RemoveType( 'TransferFunction' );
  end proc;

  zeros := proc( f:: TransferFunction) ... end proc;

  ModuleLoad();
end module:

You have to write the 'handler'; see the AddType help pages for details (it wasn't clear to me what your data structure is, otherwise I'd have suggested something). ,

You'd use seq something like this:

display([ seq( plot([X1[k],Y1[k]],[X2[k],Y2[k]]), k=1..20) ], 'insequence');

Note that your typing of v is not valid, try executing the procedure without compilation and kernelopts(assertlevel=2). Also, ?Compiler,Compile states:

Array Support
- Arrays can only appear as parameters of a procedure that is to be compiled.

While not ideal, you could make v a parameter:

kernelopts(assertlevel=2);
                                       0

v := Vector([1$4],datatype=float[8]):
foo := proc(x::float, v::Vector(4,datatype=float[8]))
    return x*v[2];
end proc:

foo(1.0,v);
                                      1.0


fooc := Compiler:-Compile(foo);
fooc := proc()
option call_external, define_external(_ma78ea510b1b4af6a145704a5af10ab6b,
MAPLE, LIB = "/var/tmp/joe-30090/_ma78ea510b1b4af6a145704a5af10ab6b.so");
    call_external(0, 1073846032, true, args)
end proc


fooc(1.0,v);
                                      1.

You can also use a regular conditional: SKIPLINES := true:

if SKIPLINES <> true then
...
end if;

to skip lines of code.

One way is to symbolically evaluate the function to generate an expression, then use unapply to create a new function:

g := unapply(f(A,'b'),'b');
                              g := b -> 3.5 + b

Note the forward quotes around b to delay its evaluation. This method won't always work, since it isn't always possible to symbolically evaluate a function. A more general technique is the following:

g := subs(_A=A, b -> f(_A,b));
                                   g := b -> f(3.5, b)

Another way to do that is the following

g := (a -> b -> f(a,b))(A);
                                   g := b -> f(3.5, b)

It is briefly mentioned in the the help page for rtable. Its first argument is a range (0..1), its second argument specifies the percentage of nonzero values (1.0, you want all to be nonzero). The elements of the rtable are randomly and uniformly initialized between 0 and 1. Using frandom is much faster than other methods (see generating an array of random floats). Here is an example of creating a Vector with 10 random elements:

rnd := rtable(1..10, frandom(0..1,1.0), subtype=Vector[column], datatype=float[8]):
                                                  [ 0.530484209985166988 ]
                                                  [ 0.952299352583890979 ]
                                                  [0.0354787515261802966 ]
                                                  [ 0.240688559954775894 ]
                                                  [ 0.767101828106371397 ]
                                         rnd :=   [ 0.724353985897650987 ]
                                                  [ 0.388089818893285866 ]
                                                  [ 0.219477841874477498 ]
                                                  [0.00149287817728582796]
                                                  [ 0.436224435569102509 ]

You can use the frandom initializer to rtable to quickly create a Vector of floats. In the following I'll generate a Matrix (rather than a Vector) so that I can use the PlotHistogram procedure from ImageTools to generate the histogram.

rnd := rtable(1..10\000, 1..1, frandom(0..1,1.0), subtype=Matrix, datatype=float[8]):
plt := ImageTools:-PlotHistogram(rnd,10):
display[plots](plt, style=point);

Just do ans mod d*e. Actually, because the inline mod operator can be reassigned to be either modp or mods (see the help page for mod), it is frequently better to use the modp or mods function directly. Thus, change the last line in the procedure to

modp({...}, d*e)

, where the ellipsis indicates the original statement (depending on which version you used.