I see that you keep bumping this Question to the top of the stack. Please don't change the Products list in the header again, unless it's to select your version of Maple. Your Question is about Maple alone; checking off other boxes is disrespectful.

I think that your problem is too general to be solved by a little bit of Maple code. Perhaps it could be the subject of a whole package with hundreds of lines of code. In particular, your iteration requires as an initialization the solution of a PDE, and it requires two indefinite integrations at each step. These are difficult things even in specific cases, let alone with the generality that you specify here. After the PDE and the indefinite integration are handled, the rest of your iteration is straightforward.

@Markiyan Hirnyk 

It's just a triviality based on the random function used, which should have no effect on the OP's problem. However, by randomizing the initial point (away from the default (0,0), which is likely a critical point), the example can be improved. I updated the code to do this, so please try again.

I have no idea where all that Typesetting junk that you show comes from. Let's see if it goes away. If not, try 1-D input.

@Kitonum I'm somewhat skeptical about using algsubs in this context, even though I don't have at hand an example where it doesn't work. I'd use simplify with side relations:

subs(t= a^4*b*c, collect(simplify(S, {a^4*b*c= t}), t));

or

thaw(collect(simplify(S, {a^4*b*c= freeze(a^4*b*c)}), t));

@emendes Please show me an executed example of what you mean by "the extra messages that come with Vector". You can attach a Maple worksheet to one of these posts. If you are trying to use op on a Vector, then that is the problem, for the same reason that nops won't work.

If you have an ordered one-dimensional structure whose contents don't need to be modified, then it's usually best to store it as a list rather a Vector.

A help file just starts as a Maple worksheet. Typically one would write one worksheet for each proc. Then use the commands in the HelpTools package to collect the worksheets, index them, and turn them into a help database.

@nm Fortunately, it's very difficult to write large Maple programs in anything other than plain text, and there's no good reason why one would want to. The majority of the Maple Programming Guide reflects this.

@emendes Note that interface(rtablesize= ...) only affects how rtables (Vectors, Matrices, Arrays) prettyprint in the GUI. It has nothing to do with programmatic access to these structures or how they print with lprint or printf.

I have interface(rtablesize= 25) set in my initialization file. You probably don't want a 1000x1000 Matrix to fully display in your GUI.

@emendes The most general command for the number of elements is numelems. It works for new- and old-style vectors, matrices, arrays, tables, lists, sets, and strings.

The operands of a new-style Vector (or Matrix, Array, or rtable) are details refering to its internal representation (such as the dimensions, datatype, indexing, storage order); nops stands for number of operands, and it returns a count of these details.

LinearAlgebra:-ColumnDimension is probably too specific, and in the long run, you may regret that you used it. It also bothers me that these low-level operations on Vectors and Matrices are not linear algebra, so they shouldn't be part of LinearAlgebra.

Why are you using row vectors? It's generally easier to use the default column vectors. It you're choocing row simply for display reasons, considering transposing before display. If v is a Vector, than its transpose is v^+ or v^%T (your choice). If LV is a list of vectors, then you can use LV^~%T.

@emendes Do you mean that the function is literally named anglebracket? Or is it one of `<,>` or `<|>`?

@emendes Given that inner is undocumented, there's no guarantee that it'll always be available. However, given its simplicity and the fact that it's used in Maple library code, I'd say that it's likely to be around. If not, it's a one-liner if you need to rewrite it. In Maple 2016:

Inner:= (u,v)-> add(u*~v):

Why do you say that it's Maple that's wrong rather than the other program? The points of intersection of the two graphs can be computed exactly and then evaluated to an arbitrary number of decimal places. I have done so (using 50 digits), and I can find no mistake in the Maple plot.

@sand15 

Use c__p/c__v. The p will never be evaluated that way. The name c__p doesn't contain the name p, just like cp doesn't contain p. The double underscore isn't an "operator". To Maple's kernel, it's just another two characters in a name. It's only purpose is that it's displayed differently by the GUI.

Regarding the syntax x__||(1..3): It's more flexible than the seq alternative: It can appear on the left side of the assignment operator, and the numbers can be letters.

@tomleslie You can change the default way that Arrays print with a `print/` procedure, like this:

`print/Array`:= ()-> `print/rtable`(<Array(args)>);

You can put that in your initialization file.

This is documented, albeit minimally, at ?print in the last paragraph of Description and the last example.

Note that you don't need to explicitly invoke a `print/` procedure. It'll be invoked automatically whenever the GUI has an Array to display. That's why I say that the procedure changes the default.

@Christian Wolinski That syntax was phased out about 20 years ago. The `.` is now the noncommutative multiplication operator, such as for matrices. The concatenation operator is now `||`. If you replace your `.` with `||`, your code will work fine.

@acer Also cat(X__, 1..3) or X__||(1..3).

@Vic ^%T is indeed transpose. But you've incorrectly transcribed what I typed. It should be ^~%T, which maps the transpose over the list that it's applied to. The reason that I used it is that I wanted the solution vectors to print horizontally. You can omit it if you're happy with them printing vertically.

Try this:

try
     SolLAM:= [LAM:-LinearSolve(2, LAM:-Mod(2, Ab, integer[]), 1, inplace= false)]^~%T
catch "matrix is singular":
     SolLAM:= NULL
end try:
SolLAM;

Only the first vector in the list is a solution. The other vectors, if any, are a basis for the null space.