You'll need to post the complete code, that is, the definition for G1, G2, and Loss.

@GPY Can you post a photo, a link to a photo, or a direct link to a webpage or PDF that has the formula from Cheney and Kincaid? Meanwhile, I'm about to post some code showing order of convergence for Halley's method. The code can be easily modified to show the order of convergence for any iterative root-finding algorithm.

The pseudocode that you give appears as nonsense to me. I can't see any relation to Muller's method. See the Wikipedia article "Muller's method".

The iterator that you give does not converge to the root at 0 for any starting value other than 0. It does converge linearly (and a very slow linearly at that) to the root at 1 for any positive starting value that I've tried. Are you sure about that iterator? I've looked through the entire Wikipedia category "Root-finding algorithms" and I can't find it. The only method that I know of with cubic convergence is Halley's method:

x[n+1]:= x[n] - 2*f(x[n])*D(f)(x[n])/(2*D(f)(x[n])^2 - f(x[n])*(D@@2)(f)(x[n]))

This has a denominator similar to yours, but there's no square root and there's the 2 coefficient. There are a few methods that use square roots, but AFAIK none has cubic convergence.

Your function has no simple roots. It only has the multiple root at 1.

@Przemek 

Like I said, for inline plots, curves will always be on top of surfaces, regardless of the order that they are listed. The disks are surfaces. It may be possible to change this order for directly exported plots, especially if you use one of the human-readable export formats such as eps or png. For jpegs, it should be possible to do something with Maple's ImageTools package, and the results can even be redisplayed inline. If you don't mind using jpeg, this is the first avenue that I'd explore; but jpegs do not make great plots. If you supply the full code, I'll work on this.

There are better solutions below. This Reply is just to show you your code debugged.

restart:
lista:= [[1,10],[2,20],[3,30]];
for i to nops(lista) do
     (x[i],y[i]):= lista[i][]
od;
(x,y):= map(convert, [x,y], list)[];

@erik10 Have you tried name quotes: `f '`[0]? Note that I place an extra space after f purely for typographical reasons.

@pimp88z Sorry, the part of my Answer mentioning Vectors was wrong. This particular command only takes data in list form. This is very unusual, especially for a command in a LinearAlgebra package, which is why I made the mistake. So, just use lists, as in

Student:-LinearAlgebra:-LeastSquaresPlot([1,2,3], [4,5,6]);

@GPY Use the restart command, then execute the above code again. Your error message indicates that at least one of y, x, or y(x) has an assigned value from some earlier code. The restart will remove that assigned value.

Note that you originally said the error message said "y(x) is expected to be of the form..."; this time you said yx. Only the latter makes sense. It is also important to report the last part of the error message, the part that begins "but received....". It was only seeing that part that led me to conclude that your error was due to inappropriately assigned values.

Your expression has no Laplace transform parameter, which is usually denoted by s.

@taro 

Thank you for inquiring about my code. For you, I am happy to explain it all.

Let M be a Maple Matrix. Then op([1,1], M) is the number of rows of M, and op([1,2], M) is the number of columns. The transpose of M is M^%T. The transpose is needed in this context because convert(M, list) lists M in column-major order, unless the Matrix has option order= C_order, which is not the default. (It seems strange that convert(Matrix(2,2, [1,2,3,4]), list) produces [1,3,2,4], but it is nonetheless true.)

In this context, $ is the repetition and sequencing operator. So, "X" $ 4 is "X", "X", "X", "X", and cat("X" $ 4) is "XXXX".

My code can be extended to an arbitrary number of columns (like Kitonum did)---and the extra comma at the end of each parenthetical grouping can be avoided---like this:

string1:= sprintf(
     ((m,n)-> cat(cat("(","%a,"$n-1,"%a)==")$m))(op(1,pt3)), convert(pt3^%T, list)[]
)[1..-3];

I find it difficult to produce code that is both brief and elegant when converting a sequence to a string. One big stumbling block is that the number of separators is one less than the number of sequenced items. You can see above that I avoided the extra comma by never including it, but that I avoided the extra "==" by including it and then removing it with [1..-3] at the very end.

Let me know if you need any more explanations.

@Preben Alsholm I get that the overall integral is divergent even if y=0, and setting y=0 is an easy way for the OP to clearly see that it is divergent.

Float(infinity) indicates that you were doing a floating-point operation such as evalf and the answer was infinity. Your integral's leading term has order 1/u^3, and int(1/u^3, u= 0..10) clearly diverges to infinity. You will easily see this if you set y=0.

@wo0olf 

Don't you think that it would be helpful for someone trying to answer your Question to see your PDE?