Your goal is not unusual, but the ways in which you've attempted to do this so far are a little unusual.

I'll show you how you can get it to work using a worksheet to build the whole module, and using the read command to get the individual command into that module. But please read the end, where I state what I consider to be a better methodology.

Your .mpl file starts out as a plaintext file. It's your canonical source for the procedure. Assemble the module by using that plaintext source. Leave out trying to use savelib or save, and don't do this:
  save(Rtest, "Rtest.mpl");

Indeed, I can think of no good scenario in which using the save command is useful or most appropriate for this task.

I'll use the same plaintext file in one of your Comments. It contains just this:

Rtest:=proc()
  print("this is a test");
end proc;

Note that when you use the read command it assigns to the global name Rtest rather than the module export instance of that name.

Here is a worksheet in which that methodology can be salvaged.

Download ronan_method1.mw

As you can see above, the process is now in two parts. The read call causes assignment to the global name, and the module definition makes an assignment to the export's name (by copying the value assigned to the global name).

Note: It's conceivable that your individual procedure had been previously stored to a .mla Maple Library Archive. I consider it far better to use the original source instead for building the module. But if you really, really want to use the instance of it from the .mla (instead of from the plaintext source file) then 1) you don't need to use the read command, and 2) you'll still need to use eval to assign the actual procedure to the module export/local.

Now here is a better way, IMO. It uses a plaintext source file for the module definition, as well as for the individual procedure(s) which form its exports or locals. It uses the $include directive rather than the read command. Here is the contents of that plaintext file:

Rtestm:=module()
  option package;
  export Rtest;

$include "Rtest.mpl"

end module;

And here is a worksheet that can read in that module definition.

Download ronan_method2.mw

I claim that this is a better way for a few reasons.

1) You didn't run into this problem with the read approach because of its simplicity, but if you have multiple module locals and exports which call each other then "fixing" the locals/exports with global name assignments can become highly problematic.

2) All your source code is now in plaintext files, which can be managed with revision control applications, and are also far more robust with regard to corruption.

3) you can use the external mint program (shipped with Maple) against the whole procedure at once. Eg, from your build worksheet, something like  ssystem("mint Rtestm.mpl"); 

So you should now be able to build your module from source. And in the same session (or worksheet if you wish) in which you've built/defined your module you can store the module in a .mla Maple Library Archive using the LibraryTools:-Create and LibraryTools:-Save commands. (IMO those are far less error prone for non-experts than is the savelib command.)

Download approximating.mw

By default that time(...) command will measure a sum of the changes of cpu time in all relevant threads. In that case you will not be able to recognize any improvement (even if it occurs), and might even misinterpret the results as denoting a slowdown.

If your machine is otherwise unloaded then you'd make easier sense from the results of the time[real](...) command, which measures wall-clock time.

Or you can wrap with the CodeTools:-Usage(...) command, and get both of those kinds of timing measurements displayed.

Did you try something like,

plots:-display(
  plots:-surfdata( Op1, xrng1, yrng1 ),
  plots:-surfdata( Op2, xrng2, yrng2 )
);

Download colorized_output.mw

Since you assigned a and b with numeric values near the beginning, don't forget to unassign them before you call NonlinearFit with an expression in which a and b are intended as unassigned parameter names.

Your candidate expression (formula as first argument passed to NonlinearFit) has a call to Dens(u/1+a*u). You get better control over what happens throughout, IMO, if you make Dens return unevaluated if its own argument is not numeric. That way you can completely seperate the working precision stuff going on within Dens from that at the higher level, without having to switch to the so-call operator form calling sequence.

Also, it doesn't seem like a good idea to force a lower working precision within the Dens procedure (with your evalf[8] call). That may confuse code within NonlinearFit that is attempting to approximate gradients at a greater working precision, and it may accidentally decide the the gradients are zero because the function isn't changing for some small change in u. (When that kind of thing happens you'd see a warning, like "no improved point could be found", or "...because first order conditions met", etc).

I ran the following attachment in 64bit Maple 14.01 for Linux.

Download nonlinearfit.mw

Note that the NonlinearFit command does local optimization, so while the computed a and b might provide a local minimum there may be a better global optima. This is one reason while the DirectSearch package (from the Application Center) is popular.

Is it possible that your data might be suitable for logarithmic fitting?

restart;

eq:=a*x^2+b*y^2+c*x*y:

terms:=[op(eq)];

                         2     2
            terms := [a x , b y , c x y]

conds1:=a>0, b<0, c>0:

conds2:=x::real, y::real, x<>0, y<>0:

select( u->is(u>0), terms ) assuming conds1, conds2;

                         2
                     [a x ]

select( u->is(u<0), terms ) assuming conds1, conds2;

                         2
                     [b y ]

remove( u->is(u<0) or is(u>0), terms ) assuming conds1, conds2;

                     [c x y]

I suggest using simple bounds (ie, ranges) for the variable, rather than adding those in as constraints. The solvers actually accept and handle those slightly differently. Using simple bounds with LSSolve I was able to avoid having to use feasibilitytolerance to guard against negativity.

Also, using tighter bounds (for all variables) allowed a better solution to be found.

Also, by using NLPSolve rather than LSSolve things went a little easier, it seemed. Note that if your simple bounds are of the symmetric form var=-a .. a then it may take 0.0 as initial point, and for your example that might result in initial gradients being zero (at which point, alas, the solver stops saying it cannot improve the initial point). I generally use a range that is not symmetric about zero gradients, if it's easy to see.

This might also be a good time to mention that the commands from the Optimization package do local optimization. If your nonlinear optimization problem has multiple local optima then it is not guaranteed that the computed solution is the global optima. See the DirectSearch package (version 2, not 1) from the Application Center, for multivariable global optimization.

I used Maple 2016.0 for this, as that seems to be your current version. I recommend upgrading to the point-release 2016.2 (which won't cost you).

Download worksheet_help_1_modif.mw

Why not use the FileTools package for opening the file? That package has commands specifically for testing whether a file exists, or is currently open, etc.

See the help page for topic FileTools,IsOpen 

Using Digits:=15 my 64bit Maple 2017.1 for Linux obtained a solution with a pretty good residual, in not so bad time.

Download bigSys_dig15.mw

Try adding the option,

tickmarks=[piticks,default]

to your call to implicitplot.

You might also wish to force the full horizontal view from 0 to Pi/2, say by also adding the option.

view=[0..Pi/2,default]

Execute the statement,

randomize():

at the start of your worksheet.

Q:=Matrix(Q, datatype=float[8]);

or,

Q:=Matrix(5,5,(i,j)->simplify(Q[i,j],zero));

or several other variants.

Similarly,

H:=Vector(H,datatype=float[8]);

or,

H:=Vector(5,(j)->simplify(H[j],zero));

or variants. As you've seen in other message threads, small nonzero imaginary components can be handled with `fnormal`, before calling `simplify`. Or you can be more forcing and use the `Re` command. 

nb.Without running your code in Maple, your Matrix K looks purely real, symmetric, and your Matrix M looks purely real, symmetric, positive-definite. If so then the eigen-solution should have all-zero imaginary components.

Are you using an older version of Maple? If I recall then in recent Maple version the zero imaginary components of complex[8] rtables are suppressed on printing.

The fnormal command let's you replace small floating-point values with 0.0 or 0.0+0.0*I, etc. This command allows you to control the fineness, optionally. See its Help page.

As Kitonum has shown, the command simplify(expr, zero) removes the 0.0 occurences (real or imaginary) from the expression `expr`.

Are you looking for plots:-implicitplot3d, where its first argument might be f(x,y,z)=K for K numeric?