@Magma You say "and also the determinant of A should be non-zero" but I'm not clear whether than means it follows automatically or it has to be another test.

I would have thought if you have an irreducible polynomial of degree n which doesn't have zero as a root (such as the one you gave), then the determinant is automatically going to be nonzero? So you only have to test if the characteristic poiynomial is equal to the required polynomial, and then you are done (or could do the determinant test only when the polynomials match to improve efficiency).

@radaar I recently used the Statistics method for a Monte Carlo method of estimating parameter errors in nonlinear regression and it seemed to work well. 

But I don't understand - your timing tests suggested it was faster than rand() and now you are saying rand is faster? Any of them can be used for Metropolis, surely.

Why don't you think the Statistics package method is random? This is what I would use.

Nice worksheet!

To simplify the log of the product to a sum you need to tell Maple N::posint (not just integer). The following works:

LL := expand(log(L)) assuming positive;

I'm assuming in the first part you didn't want to use the Statistics package, but to avoid typing errors for the PDF, I might have done

Statistics:-PDF(Normal(mu,sigma),x);

@Magma A final suggestion might be to look into the literature on totally positive matrices (all minors nonnegative) or strictly totally positive matrices (all minors positive), where there are improvements on how many determinants you need to check - for example you only need check n^2 minors for totally positive matrices, not all of them. These originate in relationships between various determinants, which might help for the MDS case though it seems the MDS case will likely be more difficult. On the other hand, all strictly totally positive matrices are MDS, so if you just want to find some MDS matrices, this might be helpful.

A good starting point is J. Pen~a, SeMA Journal, 62 (2013) 61. 

http://dx.doi.org/10.1007/s40324-013-0008-z
 

@Magma My point was mainly that doing the inverse is slower than the adjoint, and you can keep everything in integer arithmetic. And you're right that mul is probably not the best way to check all zeroes. But thinking about it, adjoint has to do the same determinants anyway, so I think the inverse idea really is not any better. I tried to do it recursively (often efficient in Maple), but it took forever.

But note that Determinant has the possibility of using many different algorithms, so you could optimize the routine to some extent. And see ?LinearAlgebra,Generic,Determinant for working in different fields. 

@Magma If Determinant(A) is 0, then it's not MDS and you are done, but otherwise just calculate the Adjoint rather than the inverse. Then for an integer matrix the Adjoint is another integer matrix. Multiplying all the entries and checking if the answer is zero is also integer-only operations. So for 16x16 this is certainly feasible. Surprisingly, mul(Adjoint(A)) is much faster than working out the Adjoint and then using mul in a seperate step [edit: no, the display just slows things down], and takes about 0.02 s on my laptop. 

 Then you can work down to the smaller matrices.

Adjoint.mw

@Thomas Richard that's embarrasing; its years since I made that mistake in my own worksheets. Perhaps pi should be prettyprinted in grey and Pi in blue, since it is a common error.

@vv taylor can do this, but I agree it shouldn't. I guess it just passes it to asympt.

Download asympt.mw

@eslamelidy plot can be done as on the attached (adjust imax to suit).

final4.mw

@eslamelidy If you use a variable b[2], then this is part of a table called b, so you cannot then use b as a variable name. I fixed this and showed how to get the table outputs you want. I'm assuming you want all these quantities at z=500.

However, you have not properly converted the integrals. The integrand of M__2 had only the integration variable t (and not z) so I used eval to substitute t for z to make the differential equation (this was just to avoid making editing errors and save me time, you could just enter the z's by hand). I then had a differential equation for M__2. But then this gets multiplied by z in creating M, which is a combination of M__2 etc. But there should be no differential equation for M. You have somehow combined into a single differential equation for M, but this has z and t together, and when you substitute t for z you just get zero, which you can see in the output of dsol3: M(z) and Mc(z) are always zero.

So you only convert M__2 etc that were integrals to differential equations, not M or Mc or N or Nc.

 final3.mw 

@eslamelidy Yes, in my method the integrals (M[1] etc) are calculated at the same time as the (other) differential equations, which makes it efficient. (You can compare these with the integrals you had after if you like.)

@eslamelidy I've started to change the worksheet, and think that if you do the rest it should work, Just to repeat, you don't need the integrals if you have the corresponding ODEs, for example f(z) = Int(sin(t),t=0..z) is the solution to diff(f(z),z)=sin(z) with initial condition f(0)=0.

final2.mw

@eslamelidy You have done a lot of cutting and pasting to get the "known" functions H(z) etc, which presumably are related simply to M[1] etc, but this is not clear in the worksheet you have given, so it is not obvious to me how to implement the details of the idea I am suggesting. You need to clearly specify how H(z) etc were constructed from the simpler functions.

I don't understand your remark.

@Carl Love I'm not seeing that in 2017.3. Executing the first time doesn't give any indication, but the second time I select and "!" (or ctrl-E) the highlighted line moves to the line with the error but no underline appears. Perhaps I'm missing a reqruired setting. I can live with the highlighted line.

CodeEdit.mw