Your model equation doesn't account for the fact that there's a vertical shift in your data. The horizontal asymptote isn't the expected 0 but rather 0.0712.... Any attempt to fit the data to the model equation will be thwarted by this.

@shzan 

Please post a worksheet showing your original code not working so that I can figure out what's going wrong.

@Bendesarts 

2*Pi isn't found because the x-values of the graph do not reach 2*Pi. The last x-value used is about 2*Pi - 1e-13. This is an anomaly of how plot chooses the x-values.

@Bendesarts 

Determing the "zero-crossing values" of an existing plot of a curve is trivial. Here's a procedure for it.

ZeroCrossings:= proc(A::Matrix)
# Input: An Nx2 Matrix extracted from a plot, representing a curve
# Output: A list of x-values of the linearly interpolated zero crossings of the curve
local Z:= table(),k;
     for k to op([1,1], A)-1 do
          if A[k,2]*A[k+1,2] <= 0 then
               Z[k]:= (A[k,1]*A[k+1,2] - A[k+1,1]*A[k,2])/(A[k+1,2]-A[k,2])
          end if
     end do;
     convert(Z, list)
end proc:

P:= plot(sin(x)+sin(2*x), x= 0..2*Pi):
A:= op([1,1], P);

ZeroCrossings(A);

     [0., 2.09456897293401, 3.14159268039291, 4.18860222947577]

@farahnaz 

Okay, I guess that's a reasonable reason to use ApproximateInt. But if you use ApproximateInt, it's up to you to adjust the parameters to control the error.

Do you have some theoretical reason to believe that the answer should change with f? If you just leave f symbolic (i.e., never assign it a value) you'll see that the final polynomial does contain f. But it contains it (apparently) in a way such that the smallest positive root doesn't depend on it.

What do you mean by "trepan"? In English, it means to drill a hole into the skull for surgical reasons. It has no other meaning, metaphorical or idiomatic.

@Bendesarts 

I think that you may be confused by Preben's imprudent use of t as the parameter. Maple doesn't care if you reuse t in this way. Perhaps you'd understand better if you replace Preben's

[Pz[i],t,t=-0.1..0.1]

with

[Pz[i],y,y=-0.1..0.1].

If you accept the easily proved fact that f(n) < < n for n > 999, it follows that h(N) = h({1..999}). Maple can easily show that h({1..999}) = {1, 2, 4}, as I'm sure you're aware.

@Rouben Rostamian  

If you're interested, here's some shortenings of your code, each procedure becoming a one-liner.

Download happy-numbers_mod.mw

@farahnaz 

My comments weren't intended to be a fix for your problem. They were intended to make tracking down the problem easier.

Why do you use evalf(ApproximateInt(...)) rather than the more normal evalf(Int(...))? The latter will return a much more accurate result. The former doesn't respect the accuracy specified by Digits.

You mentioned dsolve, but I can't find any dsolve in your code. If it's there, please highlight it with a comment.

@Kitonum RREF stands for reduced row-echelon form.

It'd be nice if you'd use comments to indicate in your lengthy worksheet where the reader can find f or D5. Also see my Reply in this other recent thread (indeed, it looks like you may be working on the same problem as that poster): "error in calculate Determinant of matrix". My point 4 in that Reply doesn't apply to you, but the others do.

Your worksheet is difficult to read and to debug for several reasons:

1. First and foremost, separate the commands into separate execution groups.
2. Use a lower value of Digits until the bugs are worked out.
3. End commands with semicolons rather than colons so that you can see their output. The infinity must occur before you take the determinant. If you can see the output of the commands, then you'll see where it occurs.
4. Do not make assignment statements of the form f(x):= y in 2D input. If you're defining a function, use f:= x-> y. On the off chance that you're making a remember table assignment, I advise switching to 1D input. It's not fair to ask the readers of your worksheet to decide between the two, and there's no way to set a default answer to the question.

If you make these changes, I'll be happy to take another look at your worksheet.

@sideshow Your attempt to use or is nonsense to Maple, although technically it doesn't violate any syntax rule, so there's no error message (until you try to use the function!). The prefix form `if` is an operator with exactly three operands. To use it for an elif situation, you need to nest them:

g:= (x,y)-> `if`(x=0.5, 200*y, `if`(y=0.5, 200*x, 'procname'(args))):

For a short chain of conditions, nested `if`s are shorter to type. But for a long chain of elifs, the parentheses required make that unbearable and you're better off using the long form if-then-elif-else-end if as in the Reply above.

But by applying logic, you can often combine conditions. For example, the four conditions from the previous Reply can be combined into a single condition like this:

g:= (x,y)-> `if`(x*y=0 or 0.5 in {x,y}, 400*x*y, 'procname'(args))):

Many complicated compound conditional statements can be made more efficient (and often less readable) by using some creative logic like this.

@sideshow 

In that case, it's better to use if-then-elif-then-else-end if:

g:= (x,y)->
     if x=0 or y=0 then 0
     elif x=0.5 then 200*y
     elif y=0.5 then 200*x
     else 'procname'(args)
     end if
:

@Rouben Rostamian  

Assuming that the OP has faithfully encoded the algorithm (which I haven't checked mostly because of annoyance with 2D-input issues and also because I don't have the book), the significant difference is that the OP uses Maple's default 10-digit software-float arithmetic and Rouben uses hardware-float arithmetic (which carries about 18 digits) for the Matrix operations. It's very likely that the book's example was executed with this same IEEE 754 arithmetic, which is called "double" or "double precision" in many languages.