@Kitonum You need to explain where you got the "magic numbers" 115, 257, and 231, because they would likely be different for other curves. I can see where 257 and 231 came from, but this may not be clear to the OP. You can always use -1 for the last index.

@Carl Love There's another sophisticated memory-management technique, which I didn't mention earlier, called weak tables. These are new to Maple 2019, and they're documented on help page ?table,weak. Basically, an item stored in a weak table gets garbage collected when there are no pointers to it other than the table itself. There are a few distinctions between items stored in the indices and those stored in the entries that you should read about on the help page.

As far as I can tell, there's no way to make the automatic remember table of a procedure a weak table, although it seems like that would be a great use for them.

Thanks, that's a clearer explanation of the problem. I think that in the other PDF, n and s were switched in some instances, and the matrix A was confused with the coefficient matrix A(t) of the IVP. Likewise, you should distinguish b from b(t) by using a differrent letter. The two vectors are not related, and not even necessarily the same size.

There's no significant benefit to making M constant rather than a function of t.

@num211 Please post this as a new Question. PDF format is fine. Any supplementary material that you can provide would be helpful (such as how to choose n). At first glance, the proposed method seems like doing a great amount of computation to get just one point of the solution (i.e., one value of t) for a linear IVP---a problem for which there are already several great algorithms, such as RKF45. 

@tomleslie I like your way of combining a colored wireframe with its corresponding surface. 

@Kitonum If your two (likely rhetorical) questions were directed at me: After making a few very minor corrections to the OP's code so that I could see the plots, I decided to color the plot's own gridlines to replace the OP's idea of drawing the gridlines with spacecurve. I think that this produced a much better plot than the OP's, while still retaining their basic structure, including coloring function. I modified the coloring function for the grid so that it would contrast. 

The OP has imaginaryunit reset to _i in their initialization file.

I like your plot.

@Kitonum Thanks. My results conform with yours if I change the initial value of t from 1 to 0. I made that initial value user-specifiable.

@Kitonum From the way that the OP has phrased the Question, it's clear (to me at least) that the solution is wanted in a matrix and vector form. They also implied (although, not clearly) that they already know how to get the dsolve solution. 

I'd have no objection to a dsolve Answer that didn't use DEtools:-matrixDE if that Answer produced a Matrix and Vector form.

@pik1432 Yes, as I said, that will work, but why do you insist on making x a matrix instead of a vector?

@pik1432 For the sixth time (!), using different words and symbols, we're telling you to make x a vector, not a matrix. The Answers from Kitonum, Mehdi Jafari, and me all show x as a vector. There is a subtle difference between a vector and a one-column matrix that is not of much mathematical significance, but which will confuse Maple if you don't respect it.

Specifically, this is what's happening: Making x a matrix will cause eq1 and eq2 to both be one-column matrices. If you index a one-column matrix using just one index, as in your eq1[i], the result is a one-element vector, not a scalar. The is command is only prepared to compare scalars, not vectors. Or, you can do the is like I showed in the Reply to my Answer below. That way will work for both vectors and matrices.

@pik1432 

S:= seq(diff~(v,W), v= V);
S[1]; #the first matrix
S[2]; #the second

@amirhadiz You must have a Maple version earlier than Maple 2018. In that case, use Kitonum's Answer instead.

@Carl Love The gradient can be also can obtained with diff~, which is similar to your original map(diff, ...). The equality of two matrix or vector expressions A and B can be verified by andmap(is@`=`, A-B, 0) (there are a great number of variations of this equality verification). 

restart:
n:= 2:
X:= Vector(n, symbol= x);
W:= Matrix(n$2, symbol= w);
G:= diff~((W.X)^+.(W.X), X); #gradient
andmap(is@`=`, G - 2*W^+.W.X, 0);

Notes:

1. Vectors and Matrices with symbolic indexed entries can be created with option symbol.
2. The transpose operator can be written ^+ rather than ^%T in some input modes.
3. The diff~ operation is equivalent to map2(diff, (W.X)^+.(W.X), X).
4. The above code requires no packages---neither LinearAlgebra nor VectorCalculus.

@awass The page ?plot,options would be too long if there were elaborate or multiple examples of each option. So, many of the options have links to supplemental help pages. The info on the linked page ?plot,typesetting applies to titles as well as to captions. This is stated in the second sentence on that page.

@mehran rajabi For example, 

eval(B, [t[0]= 0, t[1]= 1]);

This creates a new Matrix, without changing B itself; i.e., it doesn't operate inplace.

By the way, it's not clear from your original Question which t you want to change to t[1] and which to t[2] and whether you only want two new variables.