This comple tomleslie's reply

51 linear equations is not that big for Maple....
... but the output will be very large (particularly if you use w's expression and if you use integers or rationals instead of floats
So the first thing to do is not to define w! (this will be done after simplifying the solution in an ad hoc way).

Look to this attached file (no particular precaution taken to enhance the performances)

Download example.mw

1. Search Assistant in the help pages
2. Select the item worksheet, documenting, assistantsmenu (assistant)
3. Select Back-Solver
4. In the Back-Solving Assistant worksheet
   1. copy one the three (I guess) formulas that are given
   2. paste it in the empty field above the button Proceed to Back-Solver
5. Click on the button Proceed to Back-Solver, this will open a new worksheet named 'untitled..."
6. Jump into this latter, Enjoy

A suggestion: before running your computation, open the task monitor and be ready to kill the application (or the mserver process) before  it all goes wrong (in some cases using the "stop" button doesn't work).

Computing the determinant of a 3x3 matrix should be an easy task, unless each element of the matrix has a complex and lengthy expression.You seem to say that is your case? If it is so, unless some simplifications occur, the determinant will be even more complex. It could even take pages and pages to be displayed?
In such a case, what is the point to compute its formal expression if you are not able to understand it?

So another suggestion: Do the operations step by step, for instance simplify a22*a33, simplify a23*a32, form their difference, simplify (or collect/factor/normal(ize)), multiply by a11 ; do the same for the other minors
For instance

Download determinant.mw

Create an alias?

 

restart:
psi;
                              psi
alias(psi=psi(t))
                              psi

# psi can be us without explicit reference to t:
diff(psi, t)
                             d     
                            --- psi
                             dt    
dsolve(diff(psi, t)=t)
                              1  2      
                        psi = - t  + _C1
                              2         

To project a function  f(x,y) defined over some (cartesian) domain [a, b]x[c, d] you need before all to have some mapping from this plane domain onto the sphere S2.
As this is mathematically impossible you can only layer a structure on S2 in an approximate way.
 
A simple case is  [a, b]x[c, d] = [0, 2*Pi]x[-Pi, Pi] (each bounded rectangular plane domain can be mapped onto this one).
Generally this "layering" introduces singularities. To avoid them f has to be x-periodic, y-periodic, and both f(x, -Pi) and f(x, Pi) must be constant (unique values at the north and south poles);

Here is a way to "project"  a function such that f(x,y) on the unit sphere S2

projection_onto_S2.mw



 

that I found elegant but which has the drawback of becoming dramatically slow as soon as N is "slightly large".
So this is just for fun

restart:
with(PolyhedralSets):
with(ExampleSets):
N := 3:
H := Hypercube(N):
V := convert(VerticesAndRays(H)[1], Matrix)

In the command

plots[implicitplot](something, P=1..4, Q=0..3, ...);

something is aimed to be an expression which depends continuously on P and Qn (or at least which is continuous wrt P and Q but on set of null measure).
In your case this something is almost everywhere equal to 

1/(1+1.5*exp(0.01*s))

When P varies from 1 to 4 and Q varies from 0 to 3 the couple (P, Q) will take a lot of values whose only (1, 0), (1, 3), (4, 0) and (4, 3) are integers.
Thus  the definition FC:= proc(s, P::integer, Q::integer) is almost surely going to generate the error message you got.
One way coud br to write FC:= proc(s, P, Q) instead.
But in this case what are you expecting to plot?

The command 

plots[implicitplot](something, P=1..4, Q=0..3, ...);

will be equivalent to 

plots[implicitplot](1/(1+1.5*exp(0.01*s[i])), t=, P=1..4, Q=0..3, ...);

for some given value of s.
Which means you are trying to implicitplot a constant wrt P and Q.
Look for instance what happens wifh this

plots:-implicitplot(1,  P=1..4, Q=0..3)

I don not understand what you are trying do plot, but I wonder if this isn't what you would want... by any chance

Download proposal.mw

In case you would like to order the rows in such a way that only one element changes from a row to the next one you can use this 

restart
f := proc(N)
  local b := convert~(combinat:-graycode(N), binary):
  local L := [ seq(parse~(StringTools:-Explode(sprintf(cat(" %.3d"), b[k]))), k=1..2^N) ]:
  local M := convert(L, Matrix);
  return 2*~M -~ 1
end proc:

f(1), f(2), f(3)

Unless you want to use Maple for a reason of your own, and because it's the type of question asked in some magazines for the New Year, I advice you to try to solve this problem by hand.

The reasonning is here 5.mw

Here is something that is far from perfect but that can give you ideas...
proposal.mw

Tip: build the plot structure that corresponds to the icon only for non-white pixels (to avoid the icon hiding the curve).
Look to the ImageTools package to find image operations

From LinearAlgebra:-SmithFormhelp page
The SmithForm(A) function returns the Smith normal form S of a Matrix A with univariate polynomial entriesin x over a field F. Thus, the polynomials are regarded as elements of the Euclidean domain F[x].

The elements of your matrix are neither univariate (x and s) nor polynomials  but rational functions
If the variable x is a typo, proceed this way.

B := eval(A, x=s):
C := simplify(denom(B[1, 1]) *~ B):
LinearAlgebra:-SmithForm(C,s)

If the element [1,1] truly contains x, just do this

C := simplify(denom(A[1, 1]) *~ A):
LinearAlgebra:-SmithForm(C,s)

VectorCalculus[Jacobian] 
LinearAlgebra[CharacteristicPolynomial] 

How is it possible to use all of my CPU cores to calculate the eigenvalues and eigenvectors of a big matrix faster?

Write a parallel  algorithm which computes the eigenvalues.
There is a lot of articles on this topic and you could begin by this one
12_eigenvalue_8up.pdf
With Maple use the Threads package

I advice you to read the help pages for diff, dsolve and odeplot 

Download odeplot.mw

Below is an mw file that provides a detailed (though far from exhaustive) approach to your problem. However, I advise you to read the following discussion before you rush to the MAPLE file

You asked "Is there any other linear or nonlinear model that could work?"
It depends on what you mean by "that could work".

When you have data like yours that you want to "resume" by a function the first step is to wonder "What kind of résumé am I looking for?".
Do I want to set the problem in a statistical framework or not?
If it is so then a necessary first step (whilst almost always unformulated) is to model statistically the situation you are concerned with.

In your case you could for instance postulate that:

1. I have an explanatory variable named X
2. I have a dependent variable named Y
3. I postulate that  Y= F[p](X) where pdenotes a vector parameter whose "true" values pare unknown.
4. I observe Y, and possibly X, with some observational errors E_Yand E_X.
5. I postulate some relations between E_Yand E_X, for instance
   • E_Yand/or E_X  are independent (aka homoskedasticity)  of Yand X(or maybe they are not aka heteroskedasticity)
   • E_Yand E_X  are mutually independant (or not)
   • The joint distribution (general case)of E_Yand E_X  is assumed to be known.

Only after this model has been explicitely written you have the right to assess the value of p with some ad hoc method (including LinearFit, NonlinearFit, Kriging, ...)
This is how any statistical model should be build.

Now, what is this funtionF[p]?
It can be :

• an a priori model which comes from some physical considerations
• an a posteriori model infered from the visual observation of several couples (X_n, Y_n)

The postures to adopt are completely different depending on wether you have a prior model or not:

• If you have one you must not doubt it. If the data do not align with this model you should begin by verify tour mesurement processes, your assumptions of points 4 and 5 above.
  Let's take an example: you drop a marble in a column with vacuum inside and you record both the time t and the position z of the marble. You know that z=1/2*g*t²: what are you going to do if the records do not align on this parabola?
  Are you going to change the model and turn all the Physics upside-down ?
   
• Without prior knowledge of the relation Y= F[p](X)you are philosophically free to do what to do to build a model which fits the data whith any accuracy you want (in your case define  F[p](X) as a polynomial of degree 51, or use statistical models such as Kriging which can produce an exact interpolator).
  Which should restrain you do avoid doing meaningless things is a prior information you could have in the variances of E_Yand/or E_X  .
   

Without knowing what you mean by  "that could work", I have taken the liberty to make these assumptions:

• You have no prior physical model that could represent the relationship between Y and X.
  Which means that all you are looking for is a model Mwhich passes reasonably close to the observations.
• You consider that  E_X  = 0.
• You consider that the observation of Y for a given value of X differs from the what M(X) says by some random error E_Y  whose true distribution is unknown of yours.
  You just assume the variance of E_Y  is finite; that the statistical propertoes of  E_Y  do not depend on the value of  Y  itself and and that errors at different sites X and X' are mutually independent.

In this situation, which may or may not be the one you have in maind, the attached file will provide you some usefull informations (at least I hope so):
It presents at least four different models which approximate, or even interpolate your data.
Now it's up to you to clearly define your problem: what kind of model are you looking for? how do you base the goodness of a model? and so on. 

fit.mw    

All the stuff which is done in this file can be extended to ANY OTHER  type of models, even non linear ones.
I'm speaking here about the estimation of the generalization error and the selection of the "best" model


PS 1: for the CurveFitting:-Lowess solution the optimal order-and/or optimal bandwith can be asssessed with the resampling method explained in the attached file (a small bandwith will enhance the fidelity to the data and a larger one will produce a stronger smoothing but would have better predictive capabilities)

PS 2: for the CurveFitting:-BSplineCurve solution the optimal order can be found in the same way