@Rohith Like I said, you need to provide a large variety of more-complicated examples matched with the output you desire. If you don't, we're just going to go back and forth like this a hundred or so times, each time making a small adjustment to the code.

What part of the following rules don't you like? I may be able to change the rules:

1. 1 is an atomic operand.
2. v[0] is an atomic operand.
3. / is an operator.

Therefore, the only subexpression consisting of an operator acting on exactly two atomic operands is 1/v[0]. My code has precisely fulfilled the rules that you've put forth.

If you changed the input to V1/v[0] + 1, then you'd get the output that you expect. But I thought that part of the point of your project was to analyze the raw-form input, before it's simplfied (either manually or by Maple).

@taro VV has followed that help page advice 100% correctly. The problem being discussed here has nothing to do with the use of quotes.

@Kitonum The difference that this makes in the  layout of the problem is interesting. It can be simplified to

Grading:-Quiz("What is the sum?", (R,A)-> value(A)=R, ()-> `%+`((rand(1..20)$2)()), inertform);

@taro Compare the  output of these two:

eval(Diff(sin(x), x), x= 0);
subs(x= 0, Diff(sin(x), x)):

The eval knows that the x is a bound variable inside Diff, and subs does not.
 

@taro Although your subs command produces the correct result in this case, it's a bad choice for three reasons.

1. subs doesn't "know" any math, it can only blindly make substitutions, some of which may not be mathematically valid. That's why eval(e, x= a) is the preferred alternative to subs(x= a, e) in situations where something mathematical happens to e immediately afterwards, (which is most situations).

2. subs doesn't evaluate its result. For example, try subs(x= 0, sin(x)).

3. If x= a and y= b are independent substitions into e (in particular, the order that they're done doesn't matter), then subs([x= a, y= b], e) is far more efficient than subs(x= a, y= b, e).

Certainly the problem as stated cannot produce anything like that graph.

The graph looks possibly like the solution of a pair of nonlinear differential equations x'(t) = f(x,y,t), y'(t) = g(x,y,t). Could you be missing something like that?

Oh, and I agree with Mariusz. It's impolite to not give Thumbs Up or Best Answer to good/best Answers.

@acer My guess is that in order to get close to Matlab times for this, you'll need to implement a QR factorization that can be updated based on the deletion or addition of single columns from the original matrix. There are about 3 pages on this in Matrix Computations by Golub and Van Loan (let me know if you don't have access to a copy). In this particular algorithm, columns are always deleted on the left and added on the right (so the matrix's columns behave as if queued). Perhaps it is possible to exploit that special positioning of the deletions and insertions to simplify the QR updating (and thus get one-up on Matlab). I believe that QR is used in this algorithm simply because it is the factorization of dense unstructured matrices that is most amenable to being updated. 

In addition to that, you'll need to correct the OP's constant creation and destruction of rtables. If you've read the code, this is no doubt obvious to you, and I simply mention it for the benefit of other readers, especially the OP.

@kfli The embedded assignment is semantically equivalent to the non-embedded assignment, so I guess that my reasons are aesthetic, but also with a practical component. About 30 years ago was a great era of development of "software metrics": ways of measuring how long it takes (or how costly it is) to write and maintain code. (I am recalling this all from distant memory; I'll need to check, if possible, some details.) One striking thing that I read was that the overall coding time was more closely correlated with the number of lines of code than with the language, the complexity of the code, or the number of characters of code.

Regardless of whether that's actually true for most programmers, it's certainly true for me. So, I maintain a consistent and rigorous indenting style, I use a reasonable but not excessive amount of white space (my white space amount (minus the indents and right sides) is about equivalent to standard English text), and I strive to reduce the number of lines as much as possible. I put lengthy comments at the top of a procedure. I strive to get procedures to fit onto one screen. I can hold the entire procedure in my mind at once if it can fit on one screen.

@acer 

The algorithm under discussion is a method of numerically solving a system of nonlinear algebraic equations of the form F(X) = X, where X is a vector of unknowns. (Please excuse me if I'm stating something that you already know.) I don't know yet whether it's restricted to real numbers, continuous functions, Lipschitz continuous functions, etc. If one needs to solve the more-usual system G(X) = 0, I'm sure that you can see how to modify G so that it fits the format of the present algorithm.

Since it performs essentially the same job as fsolve or DirectSearch:-SolveEquations, those might be good things to compare it to.

@kfli If you want to solve df.gamma = fval by QR decomposition, you can do something like

gamma:= LinearAlgebra:-LeastSquares([((Q,R):= LinearAlgebra:-QRDecomposition(df))], fval);

with perhaps some adjustments if you're willing to store Q and R together in NAG format. Since the above will work even if df has only one column, I don't think that you should have special code for that case, at least not at this point in the development. That alone eliminates a great many lines of code.

@Ida2018 You are using the Natural units enviroment. If, instead, you use the Simple environment, then variables can "have" units, so to speak, and your original solve command should work as is, and the result will be expressed with units attached.

@Ida2018 I am putting my next Reply under my Answer, where it belongs.

@kfli 

MTM[ldivide](R, df) is equivalent to df /~ R (for scalar R and Vector df). All of Maple's elementwise operators (those that end with ~) have the same precedence, associativity, and commutativity as the corresponding operator without the ~. So df /~ R means "divide every element of df by R" (for scalar R and Vector df), which is clearly what you want.

So, the three lines can be replaced by

R:= norm(df, 2);
Q:= df /~ R;  #Normalize df.
gamma:= Q^+ . fval /~ R; #Note Maple's transpose operator.

In Maple 2018, the three lines can be (and IMO should be) replaced by the single line

gamma:= (Q:= df /~ (R:= norm(df, 2)))^+ . fval /~ R;

Please tell me whether you're using Maple 2018 so that I don't waste more time distinguishing it.

What you have labelled "Alt 1" and "Alt 3" unquestionably give the same results (modulo possibly a little round-off error, only with floating-point computation). Thus, we need to discuss whatever technique you're using to decide whether things are equal or unequal because you're doing it wrong.

@Ida2018 Sorry, I should've been more precise. I meant replace O__L by O__L*Unit(kN) in the equation. The other O__L, which is the second argument to solve, should remain as is.

@Joe Riel I think that time would be better spent by first doing obvious optimizations (obvious to you and me, that is) than in analyzing a profile or dump of a program that is obviously so raw. But first, I want to make sure that it works.