The problem with you original command is that your didn't use op. Your command was

ptslist := convert([1, 1], points, listlist);

It should be

ptslist:= convert(op([1, 1], points), listlist);

A somewhat more robust command is plottools:-getdata, because it doesn't rely on the "magic numbers" [1,1]:

convert(plottools:-getdata(points, element= "curve")[-1], listlist);

You generally can't rely on a 2d plot command to give you equally spaced data, although some will with some options (such as plot(..., adaptive= false)). 3d plots are generally evenly spaced.

If you want a 3d plot, then you need three Vectors of the same length, one for the x-coordinate data, one for the y-coordinate data, and one for the z-coordinate data. Make sure that you spell Vector with a capital V, unlike in your Question, or that you use the angle-bracket Vector constructor: < 1, 2, 3 >. If you want a 2d plot, then you just need two Vectors. Since you list two independent variables [p, sec], you are specifying a 3d plot.

To combine multiple plots, I usually prefer plots:-display over plots:-multiple, although the latter is sometimes better. Regardless of which you use, you'll need to pass it valid plotting commands. (So, note that display and multiple are both "meta" plotting commands: their primary arguments need to be plotting commands, not algebraic expressions.) The usual command for plotting equations in two variables is plots:-implicitplot. However, when it works, I prefer the command algcurves[plot_real_curve]. This command only works for some polynomial equations, but it usually gives a better plot than implicitplot when both plotting commands are passed minimal options. For example, plot_real_curves doesn't require you to specify ranges for the variables.

eqns:= [x^3-4*x=y, y^3-3*y=x]:

plots:-display(map(algcurves[plot_real_curve], (lhs-rhs)~(eqns), x, y));

From that plot we can see that all the "action" is in the ranges x= -3..3, y= -4..4. This information can be used to refine the plot using implicitplot, as done by Kitonum.

If you substitute numeric values for the variables in an equation, the result is still an equation! Were you perhaps expecting a boolean value? To get such, you need to apply at least the most trivial of Maple's boolean evaluators, evalb. In this case, (assuming that your solutions are exact algebraic numbers such as provided by solve(..., explicit)) you'll also need something like evala to simplify the algebraic numbers so that their equality can be verified.

And what do you mean by "both equations get the same result"? Even if the solution is valid, why should that mean that each equation has the same result? And it is also possible that the solution is invalid and yet the substituted equations are the same. For example, if both equations evaluate to 1=0 that would be an invalid solution, yet both equations are the same. So I see no significance to "both equations get the same result".

Now if what you actually want is to verify the accuracy of the solutions, here's something that you can do:

eqns:= {x^3-4*x=y, y^3-4*y=x}:
solns1:= solve(eqns, [x,y], explicit):
for s in solns1 do evalb~(evala(subs(s, eqns))) end do;

You will get {true} for every correct solution. This is IMO an inelegant method.

A more elegant verification method begins by converting the equations to expressions presumed to be equated to 0. But I will assume that your solution technique (your alternative to my solve) requires them to be equations. So I will convert them to expression form after generating the solutions.

eqns:= {x^3-4*x=y, y^3-4*y=x}:
solns1:= solve(eqns, [x,y], explicit):
Eqns:= (lhs-rhs)~(eqns): #Convert to expression form.
remove(s-> evala(eval(Eqns, s)) = {0}, solns1);

The returned list will contain all the invalid solutions. If they are all valid, that will be the empty list. You can use simplify instead of evala if you wish.

You have numerous syntax errors and some other problems. The specific problem that caused your error message is irrelevant. I am reluctant to list all the errors because there are so many and because some of them interact with some others in ways that make them difficult to describe individually. So first I'll post my corrected version of the code, then I'll make some comments.

B:= 10;
A:= Array(1..B, 0);

#b is in your code, but you never gave it a value.
b:= 2;
#You'll get more meaningful results by using larger primes.
p:= 2^20;

for i from 1 to B do
     p:= nextprime(p); #Not nextprime(i)
     a:= numtheory:-primroot(p);
     #Garbage collecting leads to more accurate timings.
     gc();
     A[i]:= A[i]+
         CodeTools:-Usage(
               numtheory:-mlog(b,a,p, method= indcalc),
               output= cputime, quiet, iterations= 2^12
         )
end do;

A;

Additional comments:

1. Remember to use := for assignment statements rather than simply =. 
2. Commands in packages need the package prefix, i.e., numtheory:-primroot rather than simply primroot.
3. convert(..., decimal) is meaningless in this context: The values returned by Usage are already decimal.
4. When timing a small operation, you need to use the iterations option to get meaningful results. This is because the resolution of Maple's clock is much, much larger than the resolution of your CPU's clock.
5. I don't see any point in specifying output= [cputime, output] since you're discarding the output.

To get {-2,3}, use

eval({x,y}, result);

If you use sets like that (they're sets because they're in curly braces), then you cannot rely on the first element corresponding to x and the second corresponding to y. Indeed, if the values of x and y are the same, then your set will only have one element. For these reasons, you should use lists instead of sets, like this:

eval([x,y], result);

The $ (unlike seq) can work symbolically, i.e., for symbolic values of i. And diff can work with this information to produce derivatives of arbitrary order. (Of course, diff can only do this for somewhat simple functions.) Observe the output of your first command:

diff(x^7, x$i);

This expression gives a valid derivative for any i, including i = 0. In your first example, this is the expression that is passed to seq, which then simply plugs in the values of i.

The command seq has "special evaluation rules" similar to the first parameter being declared with uneval (see ?uneval) . (These rules also apply to add and mul.) You are correct that these rules are not described at ?seq. These rules mean that in your second example, the unevaluated expression diff(x^7, x$i) is passed to seq. Evaluating this expression at i=0 causes an error because the expression becomes diff(x^7, NULL), i.e., nothing is after the comma. This can be avoided by using diff(x^7, [x$i]), which is valid for any i, even if the [ ] are empty.

You have a classic off-by-one error. In addition to the problem pointed out by sazroberson of not propagating in the other arrays, you are also not propagating to the next value of i.

Change the lines

S[i]:= S[i-1], etc. to

S[i+1]:= S[i], etc.

In the long if - elif block, change each line of the form

S[i]:= S[i]+1 to

S[i+1]:= S[i]+1.

Remove the if statement that prevents the execution of this code for i=0.

Jakub,

I have finished the project including extensive research of your course handouts from your professor's website (http://www.kosiorowski.edu.pl/wp-content/uploads/2014/10/Ia-NumMet-W-S-1-Error-theory.pdf) and a complete-from-scratch rewrite of your code including individually checking for zero denominators for each method (so that the procedure will perform all applicable methods even if some are inapplicable). I extended the code so that it works for functions of any number of variables (not just x, y, z) while at the same time vastly simplifying it.

Download INVERROR.mw

The expression seq(a[i_j], i= 1..3) makes no sense. I don't just mean that it makes no sense to Maple; I mean that it makes no sense at all. For one thing, you have no iteration over j. And even if you did iterate over j, you would then have subscripts on the integers i. And what would that mean?

Something that could possibly make sense would be

seq(a[i[j]], j= 1..3)

or

seq(a[i||j], j= 1..3)

or

seq(a[i_||j], j= 1..3).

In each case, the iteration is over j. Perhaps that is what you meant.

Your problem is caused by you typing an extra space after cos. You have entered the forcing function as

u*cos (omega[g]*tau)

which is interpretted by the odious 2d input parser as if you had typed

u*cos*omega[g]*tau.

You need to change it to

u*cos(omega[g]*tau).

It's only a difference of one space. 

This problem will be avoided if you use the 1d input (which ignores extra spaces between a name and a punctuation mark, like most reasonable computer languages), like this:

(Display of uploaded worksheets is broken. But you can download my worksheet below.)

Download ODE.mw

I recommend that you do not use .jpg for mathematical plots: it makes lines and text blurry. So, it'll be .png.

The plotsetup command can be used to programmatically redirect plot output to a file. I find this easier and more precise than using mouse controls to export a high-quality plot (although I do use the mouse controls to export most of the regular-quality plots that I post here). Working with high-resolution plots in a worksheet (either Standard or Classic) will make the interface excruciatingly slow. Of course, you'll probably want to look at the plot in the worksheet while you're perfecting it. While you're doing that, avoid scrolling the worksheet, and never have two such plots in the same worksheet.

plotsetup(png, plotoutput= "myplot.png", plotoptions= "width=1920,height=1080");

After issuing the above command, all plotting-command output will go to the file. To turn this off, use plotsetup(default). Notice that I used high-resolution specifications for the width and height. There is also a tranparency option that can be set to a number between 0 and 1 inclusive.

There's little point in export at high resolution if you haven't computed enough points in your plot. So, use the plotting-command options such as numpoints and grid to raise the number of computed points.

The reason that it's taking so long is that your inner integrals are symbolic, and Maple is not finishing that symbolic integration. Maple can do your innermost integral symbolically in a few seconds. The result is about three screens long. The next layer out takes longer; it may be the one that doesn't finish.

For the vast majority of Maple commands (procedures), the arguments are evaluated before being passed. For example, if your call were simplified to

int(int(f(x,y), x= 0..g(y)), y= 0..b, numeric);

that inner integral would be performed symbolically before the numeric integration is even started.

There's no point in trying numeric integration when you have symbolic parameters. You have k, l, a, and b.

Maple has an alternate simplified notation for multiple integrals that will help you avoid the problem of premature evaluation of inner integrals, and also vastly simplify the parentheses needed for your four-deep integral. This notation works for both numeric and symbolic integrals. In the alternate notation, my integral above becomes

int(f(x,y), [x= 0..g(y), y= 0..b], numeric);

Because Vector is an active procedure as well as a type and metatype, the procedure call Vector(Element) needs to be prevented. If you want to do this and also use TypeTools:-Exists, then you need two pairs of unevaluation quotes on Vector:

TypeTools:-AddType(ExpandedLine, ''Vector''('Element'));

(That's not double quotes; it's two sets of single quotes.) I can't tell you a concrete rule by which I know to use two sets of quotes; however, I've tested it, and it works.

Unlike arrays, two lists L1 and L2 can be checked for equality as a whole, i.e., without needing to do it elementwise:

if L1=L2 then "They're equal." else "They're not equal." end if;