The benefit of using a ".m" file is that it automatically stores all the information needed to reconstruct all the objects stored in it, whether they be arrays or anything else; and the entire reconstruction is handled quickly by a simple read command. But storing the array in a binary file is even much faster than a ".m" file (by a factor of 10 or so).

In the code below, my procedure SaveArray stores three files: the first a pure binary file containing the array's data, the second a small ".m" file containing the array's attributes (you don't need to know about attributes to understand the rest of this), and the third a small file of plaintext Maple code that can be used to exactly reconstruct the array from the first two files in a new session so that it's structured exactly the way that it was in the originating session.

restart
:
(*---------------------------------------------------------------------
Procedure to save a hardware-datatype rtable in a binary file along
with all the info needed for its reconstruction in a new session.

The 1st argument is the rtable, which must've been assigned to a 
name. The 2nd argument is a file name "stem" to use for three files.
It's assumed that currentdir() has been set to the desired directory.

The output is 3 files and a snippet of text on screen. For example,
if the stem is "mydata", then the output files are
    mydata.bin    | the raw form of the rtable
    mydataAttr.m  | the rtable's 'attributes'
    mydata.mpl    | plaintext Maple code file to be 'read' to 
                  | reconstruct the rtable from the 1st 2 files
The snippet of text is a line of code that can be copy & pasted to a
new session to further automate the reconstruction.
---------------------------------------------------------------------*)
SaveArray:= proc(A::evaln(rtable), F::string)
description 
    "Save a hardware-datatype rtable in binary file along with all "
    "info needed for its reconstruction in a new session           "
;
option `Author: Carl Love <carl.j.love@gmail.com> 2021-May-27`;
local 
    O:= [rtable_options](eval(A)), 
    OE:= select(type, O, `=`)
;
    (:-_Attr, O):= selectremove(
        type, O, identical(attributes)=anything
    );
    :-_Attr:= :-_Attr[];
    save :-_Attr, cat(F, "Attr.m");
    FileTools:-Binary:-WriteFile(cat(F, ".bin"), eval(A));
    FileTools:-Text:-WriteFile(
        cat(F, ".mpl"),
        sprintf(
            "read \"%sAttr.m\":                                  \n"
            "%a:= ArrayTools:-Alias(                             \n"
            "    FileTools:-Binary:-ReadFile(                    \n"
            "        \"%s.bin\", ':-datatype'= ':-%a'            \n"
            "    ),                                              \n"
            "    %a, ':-%a', %a[]                                \n"
            "):                                                  \n"
            "rtable_options(                                     \n"
            "    %a, %a[],                                       \n"
            "    attributes= [                                   \n"
            "        eval(                                       \n"
            "            attributes,                             \n"
            "            select(type, [rtable_options](%a), `=`) \n"
            "        )[],                                        \n"
            "        rhs(:-_Attr)[]                              \n"
            "    ]                                               \n"
            "):                                                  \n"
            ":-_Attr:= ':-_Attr':                                \n",
            F,                     
            A, 
            F, eval(':-datatype', OE),
            [rtable_dims](eval(A)), eval(':-order', OE),
            `if`(':-readonly' in O, [':-readonly'= true], []),
            A, remove(type, O, identical(datatype, storage)=anything),
            A
        )
    );
    :-_Attr:= ':-_Attr';
    printf("#Copy & paste next line to new session and execute:\n");
    (printf@cat)(
        "currentdir(\"",
        StringTools:-SubstituteAll(currentdir(), "\\", "/"), 
        "\"): ",
        "read \"", F, ".mpl\":\n"
    );
    return
end proc
: 
#random array simulating your conditions:
B:= rtable(
    (1..434000, 1..28), frandom(0..1), 
    datatype= hfloat, order= C_order, readonly, subtype= Matrix,
    attributes= ["This is the good data"]
):
B[42, 23]; #some random entry for verification
                       0.528682555121851

currentdir("/Users/carlj/desktop"):
st:= time[real]():
   SaveArray(B, "mydata");
time = (time[real]()-st)*Unit(second);
#Copy & paste next line to new session and execute:
currentdir("C:/Users/carlj/desktop"): read "mydata.mpl":

                      time = 0.126 Unit(s)

restart:
st:= time[real]():
    currentdir("C:/Users/carlj/desktop"): read "mydata.mpl":
time = (time[real]() - st)*Unit(second);
                      time = 0.164 Unit(s)

B[42, 23];
                       0.528682555121851

rtable_options(B);
datatype = float[8], subtype = Matrix, storage = rectangular, 
  order = C_order, readonly, attributes = [source_rtable = 
  [...numeric data elided...],
  "This is the good data"]

@2cUniverse Both the original procedure and VV's version require zero padding of the lists if either are shorter than d. I've made that correction and also added some other error checking:

`&*`:= proc(A::list, B::list, b::And(posint, Not(1)):= 10, d::posint:= 4) 
local i, k, c:= 0, A0:= table(sparse, A), B0:= table(sparse, B);
    if [A,B]::list(list(integer[0..b-1])) then 
        [seq](irem(c + add(A0[i]*B0[k+1-i], i= 1..k), b, 'c'), k= 1..d)
    else
        error "expected lists of base-%1 digits, received %2", b, [A,B]
    fi
end proc:

@Carl Love Actually, the same workaround can be done without overload because the eval pre- and post-wrappers do nothing in the usual case that the input doesn't contain a free _Z. I'd rather not use overload if it can be avoided because the structure that it returns is a bit awkward to deconstruct. Thus, 

restart:

interface(warnlevel=4):
kernelopts('assertlevel'=2):

CI__orig:= eval(SolveTools:-CancelInverses):
unprotect(SolveTools:-CancelInverses):
SolveTools:-CancelInverses:= e->
local _A;
    eval(CI__orig(eval(e, _Z= _A), _rest), _A= _Z)
:
protect(SolveTools:-CancelInverses, CI__orig):

expr:= exp(x)*sin(y)-3*x^2+(exp(x)*cos(y)+1/3/y^(2/3))*Z = 0;
solve(expr, y);

@tomleslie Ah, I was led astray by an incomplete 3rd-party help file https://exceljet.net/excel-functions/excel-mround-function which doesn't specify the type of m, and shows numerous examples using integer m. 

@MapleEnthusiast Your array has numIters rows. Row 1 stores the variable labels. Row 2 stores the initial guess. That leaves numIters-2 rows available for storing iterations. You need to create the array with 1 more row.

@mmcdara You said:

•  it is not recommended to work with csv files this size.

That depends on who's doing the recommending. If computational efficiency is the concern, then I'd agree that you don't want to be constantly rereading a huge CSV. On the other hand, if data archiving is a concern, then I'd strongly recommend keeping it as CSV, which is essentially a universal format which'll remain easily readable for many years.

To the OP: If you use either acer's or my solution to this problem, then this is not an issue. So I recommend that you don't change that file!

@kiraMou You're still confusing singular "solution" with plural "solutions". Consider the solution list4 from your posted worksheet. It happens that this solution has imaginary parts -Pi/2 for all four variables. But you think that half of them should be +Pi/2. That's not what's meant by "the solutions (plural!) are symmetric with respect to the imaginary axis." That statement doesn't mean that each individual solution (each of which contains 4 variables) is symmetric. Instead, it means that if you take the conjugates of all 4 values, then you'll get another solution. Let's do that:

list4c:= evalindets(list4, complex(numeric), conjugate);
eval({t[1],t[2],t[3],t[4]}, list4c);

Now you can see that the 4 imaginary parts have been changed to +Pi/2. And you can see that the residual imaginary parts from this new solution are effectively 0, the deviations being due simply to roundoff error, and the residual real parts are effectively -1. Thus list4c is also a solution. It is the symmetric counterpart of list4.

@Carl Love Furthermore, if one did want to represent a curve as a surface (and there are good reasons that one might want to do that, such as elaborate colorings), then the better command for that is tubeplot, as shown by mmcdara.

@kiraMou A correct assignment of numeric values to x[1], ..., x[4] is a single solution. You've been referring to such a thing as "solutions" (plural). It's possible that this misclassification has led to some of your confusion.

@nm You wrote:

• Error, (in my_module_name:-dsolve) assertion failed in assignment to ode_obj, expected ode_type, got _m1657072644000 

Okay, this is some progress, because this error message indicates (because of the strange name _m1657....) that an object module rather than just a symbol was being passed between procedures. The problem is now that it doesn't recognize that object as being type ode_type. Both immediately before and immediately after you receive this error, at the top (global) level, give the command TypeTools:-Exists(ode_type).

• Right now, everything is working fine, once I remove the ::ode_type from the local variable definition and keep the assert there also.

That's an unsatisfactory solution, as you likely realize. But the more-important thing is is that returned object of the correct type, even though perhaps that type can't be properly named?

I look forward to your MWE that truly encapsulates the issue. This problem is something that I've struggled with myself: An interconnected family of modules, some of which define types in their ModuleLoad, others with module locals declared to be of those types, running under assertlevel=2. There is some deep-seated bug in this situation, even when none of the modules are objects. 

@Gabriel samaila There is nothing in the code posted in my Answer above that won't work in Maple 16. If your original code works, then that code will work. Could you post a new executed worksheet showing it not working?

@Kitonum If you use plot3d for a curve, as you just showed, it does the computations as if it were doing a surface, producing a 49x49x3 MESH that is constant along one dimension. The black that you see is the grid lines of that mesh "mushed" together under the default style= surface. Thus, if you changed that to style= wireframe, it would respect your color setting.

Thus, I fully stand by my statement "plot3d is only for plotting surfaces", because the data structure produced is for practical computational purposes a surface, even if a topologist would formally call it a curve.

@Kitonum I wonder why I needed to add some transparency to the surface and you did not. I'm using Maple 2021.

@nm Yes, I was writing my Answer at the same time as you were writing your Edit.

You wrote:

1. The other proc called is returning back an object of ode_type....
2. things work for other types, like integer, string, etc....

Neither of those statements are true. Perhaps you should post an example that leads you to believe that 2 is true.

A declaration of the form local A::typeA; does not create an entity A such that type(A, typeA) is true. Rather it checks (under assertlevel=2) that that is true for every := assignment made to A (and it may check some other forms of assignment also). A is just type symbol until you make an assignment to it. In the case at hand, you've never made an assignment to ODE.

You'll need to be more precise about what you mean by "touch", because a circle can be found for any 3 noncollinear points. Touching could mean a shared tangent line, a shared radius of curvature, or perhaps something else.