Like I said to your earlier version of the question: Do you really want to use RK4? Is your goal to solve the ODE, or is your goal to complete a programming exercise for RK4? You will need to specify a numeric value for k_1 either way.

@Kitonum Yes, ilog10, ilog2, and length are much faster than ceil@log[10] for exact aritnmetic.

@Joe Riel The procedure has option remember, but it is not an essential part of the computation like it is with the recursive procedure.

@Andriy In your animate command, you use u(t,x). In the vast majority of cases, arguments are evaluated before they are passed.

Do you really want to use RK4? There are easier ways. To use RK4, we need to convert to a first-order system. And however it's solved, we'll need another initial condition.

@niller123 Sorry, there's currently a bug with uploading files. I'll send you an email, and you reply to it with the file.

@Axel Vogt The procedure produced by rsolve with option makeproc uses a "sliding window" such that it saves the previous three computed values. Aside from error checking, it's equivalent to this:

X:= proc(n)
     local k,L,r;
     L:= [0,1,2]; #Initial values
     for k from 3 to n do
          r:= `+`((k*L)[]);
          L:= [L[2], L[3], r]
     end do;
     L[4+n-k]
end proc:

@mois Can you show the exact command that you used and the results that it returned? Also, why didn't you do the model function f:= a+b*x^c+d*x^e?

Yes, please "send" the document (worksheet) by attaching it to a Reply.

By the way, your English is fine by MaplePrimes standards. Try to make a better title though.

As a moderator, I changed this Question's title from "A Rounding bug in maple" to "Is this a Rounding bug in maple?"

@Sujaan Kunalan There is no way to explicitly mark a Question as answered, but you can "vote up" for the Answer by clicking on the "thumb up" icon in the upper right corner of the Answer.

@Mac Dude The only thing that makes the difference in your call to NonlinearFit is the initialvalues. All Statistics package routines that take a Vector will take a list also.

@Markiyan Hirnyk Quoting from ?combine:

For many functions, the transformations applied by combine are the inverse of the transformations that are applied by expand. For example, consider the well-known identity
           sin(a + b) = sin(a) cos(b) + cos(a) sin(b)
The combine function applies the identity from right to left, whereas the expand function does the reverse.

You wrote: What is the reasoning behind the output?

The reasoning is that the result should also be valid for n=0.

@Alex Smith Why [-90,1] and not [-90,0]? Also, you can hide the z-axis label.

plot3d(
     [x,y,0], x= 0..1, y= 0..1-x,
     axes= boxed, labels= [x,y,``], style= patchnogrid, orientation= [-90,0]
);