I see three possibilities

1. Mapleprimes does not accept uploaded files
2. Your browser settings do not allow file uploads
3. Your firewall/security processes do not allow file uploads

I'll deal with each of these in turn

1. All(?) Mapleprimes users, including me, are able to upload files with no problem, so I think we have to accept that there isn't a real problem with the Mapleprimes site
2. Depending on you browser (and its settings) certain things may be "forbidden". For example, my default browser is Firefox, and it *really* does not like me copying hypertext links around the Mapleprimes site (or any other site for that matter).. I know how to switch off this Firefox restriction, but I choose not to. I can usually work around it. In order to check whether you have a browser restriction, I can only suggest that you experiment with a different browser. When I decide that I really *really* need a hypertext link in a post on Mapleprimes, then I just use either Chrome or IE11 (either of which will work because their default security options are a bit more "relaxed". So does the same problem occur with all browsers - if not then you need to adjust the security settings on your default browser (or use a browser which works)
3. If no browsers allow file uploads to Mapleprimes, then I would recommend that you try file uploads to other "similarish" websites. For example download a matlab .m file from somewhere and try uploading to http://www.use-matlab-online.com/matlab-user-forum.html#/ (you will probably have to register first). Or do the same with Mathematica. Basically can you upload any files to anywhere? If you can't then the problem is with your site security processes.

For the first part of your question, the easiest way is probably to generate two plots and then display them together, as in the 'toy' example below

ymax:=10:ymin:=5:fn:=x*sin(y):
a:=plot3d( fn, x=-5..5, y=-ymax..-ymin, color=red):
b:=plot3d( fn, x=-5..5, y=ymin..ymax, color=blue):
plots[display]([a,b]);

I don't really understand the second part of your question (and some of it is missing anyway)

Carl is correct, but actually you don't even need a display() command, because plot3d() will accept a list of plots as its first argument, so

plot3d( [seq(srr[j], j=1..n)], +whatEverPlotOptionsYouWant)

ought to work. I have tested this with

srr[1]:=r^2*sin(theta):
srr[2]:=r+theta^2:
plot3d( [ seq(srr[j], j=1..2)], r=0..1, theta=0..2*Pi, color=[red,blue]);

and it works perfectly

Assuming 'A' is some kind of sequence, then [A] will produce a list

Without affecting any assignments (but with less info), try

t1 := time();
for i to 3 do
    c[2] := Vc[2][i];
    dsol[i] := dsolve(dsys, numeric, continuation = lambda);
    print(c[2]);
    print(dsol[i](0))
end do;
t2 := time()-t1

On Maple 2016.1 go to

Tools->Options->Interface

and click the "Confirmation Dialogs" button. Then check the status of the "Save Large Calculations". You want this box checked: that *should* give you a confirmation dialog when saving a lot of stuff

You posted an (almost illegible) picture, rather than using the big green up-arrow to post code. This means that anyone responding has to retype the contents of your 'picture'. Believe me, very few people here will retype your code, just because you are too lazy/incompetent to post it in a way which allows us to use copy/paste.

Now in legible form you have

pdesys:={ diff(J(x,t),x)=epsilon*J(x,t)*C(x,t),
                diff(C(x,t),t)=phi*epsilon*J(x,t)*C(x,t)
              };

and Maple return the answers

{C(x, t) = 0, J(x, t) = _F1(t)}, {C(x, t) = _F1(x)*_C1, J(x, t) = 0}

Well let us consider these, the first answer has C(x, t) = 0, so your differential system becomes

diff(J(x,t),x)=0
diff(C(x,t),t)=0

I have news for you, if C(x,t) is identically zero, then so also is diff(C(x,t),t), so the solution is valid for the second equation. It is also true that provided J() is any function which does not depend on 'x', then the first equation will be satisfied, so specifying J(x,t)=_F1(t) (ie independent of 'x') is completely correct. It just means that J() can be any arbitrary function of 't'

Now for the second answer: if J(x, t) is identically zero, then once again your PDE system reduces to

diff(J(x,t),x)=0
diff(C(x,t),t)=0

Since J(x,t) is identically zero, then so also is diff(J(x,t),x), so the first equation is satisfied. The second equation will be satisfied provided that C(x,t) is a function of 'x' only, so that diff(C(x,t),t)=0. Hence specifying the solution for C() as _F1(x)*_C1 is completely correct (although I admit I am not sure why Maple with bothered the arbitrary constant _C1)

So the solutions provided by Maple are completely correct: are other solutions possible - maybe. ~You claim to have analytic solutions for these PDEs where neither J(x,t) nor C(x,t) are zero - please post them.

If I convert your code to 1-D math input, then I can see a few issues

1. As Carl has already pointed out you should be using pdsolve() not dsolve()
2. I think(?) the syntax of many of your boundary conditions is incorrect - see amendments below
3. Several of your boundary conditions depend on the variable 'h', which is nowhere defined

If I correct your code (and insert h:=1, just to give this necessary variable some sort of value), then I end up with the code

restart;
d[11] := 1; mu[11] := 1; q[311] := 0; d[33] := 1;
mu[33] := 1; a[11] := 1; e[311] := 0; a[33] := 1;
A := 1; g[111111] := 1; c[1111] := 1; g[113113] := 1;
f[3113] := 1; beta[11] := 1; `ΔT` := 1;
II := 1; L := 1; h := 1;
J := d[11]*(diff(Phi(x, z), x, x))
     +
     mu[11]*(diff(psi(x, z), x, x))
     +
     q[311]*(diff(w(x), x, x))
     +
     d[33]*(diff(Phi(x, z), z, z))
     +
     mu[33]*(diff(psi(x, z), z, z));
B := a[11]*(diff(Phi(x, z), x, x))
     +
     d[11]*(diff(psi(x, z), x, x))
     +
     e[311]*(diff(w(x), x, x))
     +
     a[33]*(diff(Phi(x, z), z, z))
     +
     d[33]*(diff(psi(x, z), z, z));
dsys := {B, J};
BCS := { Phi(0, z) = 0, Phi(L, z) = 0, Phi(x, -(1/2)*h) = 0,
         Phi(x, (1/2)*h) = 0, psi(0, z) = 0, psi(L, z) = 0,
         psi(x, -(1/2)*h) = 0, psi(x, (1/2)*h) = 0
       };
dsol5 := pdsolve(dsys, BCS, numeric);

However, the problem now is that they values you have assigned to the coefficients

d[11] := 1; mu[11] := 1; q[311] := 0; d[33] := 1;
mu[33] := 1; a[11] := 1; e[311] := 0; a[33] := 1;
A := 1; g[111111] := 1; c[1111] := 1; g[113113] := 1;
f[3113] := 1; beta[11] := 1; `ΔT` := 1;
II := 1; L := 1; h := 1;

mean that the differetnial equations J and B are identical, so actually you have only one differential equation, containing two dependent variables. Maple then (quite correctly) complains that

Error, (in pdsolve/numeric/process_PDEs) number of dependent variables and number of PDE must be the same

So you cannot reasonably expect to solve this system

BTW It seems rather pointless to include the dependent variable w(x), and then set a coefficent equal to zero to ensure that it 'disappears' from your equations. However if the coefficient of diff(w(x), x, x) is ever non-zero, then yuo will have to produce another pde (three dependent variables requires three equations), as well as extra boundary conditions for w(x). Since it only has a single independent variable and is of degree 2, then another two boundary conditions should suffice

Nothing wrong with Markiyan's solution, but alternatives are always worth considering.

Having produced a histogram, one can use the plottools:-getdata() command to get the underlying structure of what has been plotted. With this information it is relatively easy to produce sequences of [binNumber, binContents] or [binMidPoint, binContents], or almost anything else you might want

Consider for example the attached worksheet

bins.mw

I tidied up your code to avoid the use of 'deprecated' constructs, since these were making debug efforts more difficult. Modified code is attached here

isom.mw

This code generates the following set of equations for the coefficients of the C-matrix

0 = 0,
0 = c[2,2]^2,
0 = c[1,2]*c[1,1],
0 = c[1,2]*c[2,2],
0 = c[2,2]*c[1,1],
0 = c[2,2]*c[2,1],
c[1,1] = 0,
c[1,2] = c[1,2]^2,
c[2,1] = c[1,2]*c[2,1],
c[2,2] = 0

Given c[1,1] = 0, and c[2,2] = 0, the above set of equations can be reduced to

c[1,2] = c[1,2]^2,
c[2,1] = c[1,2]*c[2,1]

Rather obviously, this reduced set has two solutions:

1. c[1, 2]=0 and c[2,1]=0 (in other words the matrix C is identically zero)
2. c[1,2]=1, in which case c[2,1] is arbitrary, which is why you get the solution [[0, 1], [c[2,1],0]

PREDS and ZPLOT are local to the procedure Z_TABULATE( ZMAT ). Their contents are returned when the procedure is executed.  But these contents are are not assigned to anything when returned by the procedure

You can either change the procedure call to

p1, p2:=Z_TABULATE( ZMAT ) ;

so that  the contents of PREDS, ZPLOT are assigned to p1, p2 respectively, and then use

pointplot(  p1, p2 , connect , color = "Green");

(Or you could shift the pointplot() command inside the procedure)

Commands within the DynamicSystems package will allow you to examine/manipulate poles/zeros and characteristic polynomials etc of any kind of linear systems. As such you can implement almost any kind of stability test which you want

For Nyquist stability, check out

DynamicSystems[GainMargin]()
DynamicSystems[PhaseMargin]()

For Lyapunov stability, check out

DynamicSystems[Covariance]()
DynamicSystems[Grammians]()

For Routh stability, check out

DynamicSystems[RouthTable]()

Not too familiar with the Jury criterion, but it seems to depend on relatively simple manipulations of the coefficients of the system's characteristic polynomial, so

coeffs( DynamicSystems[CharacteristicPolynomial]() )

ought to be useful

you could try the basic secant method (or any of its higher order variants), as in

h:=1e-05;
diffu:= (xx, tt) -> (rhs(pds:-value()(xx, tt+h)[3])-rhs(pds:-value()(xx, tt)[3]))/h;
diffu(0.5, 0.5);

I'm not exactly sure what your problem is.

In the attached, I have added some execution groups to your worksheet which illustrate the use of the 'plot',  'plot3d', and 'value' methods

U(R)_addons.mw

By default, the independent variable, is returned or plotted. It is also possible to calculate or plot any function of the independent and dependent variables, provided that this function can be numerically evaluated at a point.

You indicate that you want (eventually?) to use d_theta(x,R)/d_R. Although this is a function of the dependent variable, it cannot be evaluated from numeric data. I suspect that the best you will be able to achieve, is to compute a sequence of values for theta(x,R) and then use some kind of successive differences to approximate this derivative

See the attached for some fixes

1. I specified a midpoint method for the dsolve() command
2. dsolve() was then unable to achieve default accuracy (abserr=1.0e-05) with the default 'mesh', so I increased the maxmesh value  (One could alternatively, keep the same mesh and slacken the error tolerancwe)
3. I added a plot() command - just to check that some sort of sensible output was obtained

1MASS_Corr.mw