Actually, we can simplify it somewhat.

> normal(fnormal(%));

{v[ds](t) = .4958026528*exp((30764029.75+1577330139.*I)*t)*`#msub(mi("V"),mo("in",fontweight = "bold"))`+9.149694069*exp((30764029.75+1577330139.*I)*t)-(.766755745*I)*exp((30764029.75+1577330139.*I)*t)-(0.3290893603e-1*I)*exp((30764029.75+1577330139.*I)*t)*`#msub(mi("V"),mo("in",fontweight = "bold"))`-0.3002643e-3*exp(-10002901.82*t)*`#msub(mi("V"),mo("in",fontweight = "bold"))`-42.51993430*exp(-10002901.82*t)+0.1702286306e-1*exp(-0.1544864054e11*t)*`#msub(mi("V"),mo("in",fontweight = "bold"))`+.43114634*exp(-0.1544864054e11*t)+24., v[gs](t) = -0.3674479877e-3*exp((30764029.75+1577330139.*I)*t)*`#msub(mi("V"),mo("in",fontweight = "bold"))`-0.8037470713e-2*exp((30764029.75+1577330139.*I)*t)-(0.7193539082e-1*I)*exp((30764029.75+1577330139.*I)*t)-(0.3899911367e-2*I)*exp((30764029.75+1577330139.*I)*t)*`#msub(mi("V"),mo("in",fontweight = "bold"))`-0.6005358003e-4*exp(-10002901.82*t)*`#msub(mi("V"),mo("in",fontweight = "bold"))`-8.504088823*exp(-10002901.82*t)+0.1556340602e-2*exp(-0.1544864054e11*t)*`#msub(mi("V"),mo("in",fontweight = "bold"))`+0.3941819611e-1*exp(-0.1544864054e11*t)+12.}

The <maple> tag doesn't recognize the typesetting codes for V_in: make that

{v[ds](t) = .4958026528*exp((30764029.75+1577330139.*I)*t)*V[`in`]+9.149694069*exp((30764029.75+1577330139.*I)*t)-(.766755745*I)*exp((30764029.75+1577330139.*I)*t)-(0.3290893603e-1*I)*exp((30764029.75+1577330139.*I)*t)*V[`in`]-0.3002643e-3*exp(-10002901.82*t)*V[`in`]-42.51993430*exp(-10002901.82*t)+0.1702286306e-1*exp(-0.1544864054e11*t)*V[`in`]+.43114634*exp(-0.1544864054e11*t)+24., v[gs](t) = -0.3674479877e-3*exp((30764029.75+1577330139.*I)*t)*V[`in`]-0.8037470713e-2*exp((30764029.75+1577330139.*I)*t)-(0.7193539082e-1*I)*exp((30764029.75+1577330139.*I)*t)-(0.3899911367e-2*I)*exp((30764029.75+1577330139.*I)*t)*V[`in`]-0.6005358003e-4*exp(-10002901.82*t)*V[`in`]-8.504088823*exp(-10002901.82*t)+0.1556340602e-2*exp(-0.1544864054e11*t)*V[`in`]+0.3941819611e-1*exp(-0.1544864054e11*t)+12.}

You could also use evalc to express the complex exponentials in terms of real exponentials times sines and cosines.

Actually, we can simplify it somewhat.

> normal(fnormal(%));

{v[ds](t) = .4958026528*exp((30764029.75+1577330139.*I)*t)*`#msub(mi("V"),mo("in",fontweight = "bold"))`+9.149694069*exp((30764029.75+1577330139.*I)*t)-(.766755745*I)*exp((30764029.75+1577330139.*I)*t)-(0.3290893603e-1*I)*exp((30764029.75+1577330139.*I)*t)*`#msub(mi("V"),mo("in",fontweight = "bold"))`-0.3002643e-3*exp(-10002901.82*t)*`#msub(mi("V"),mo("in",fontweight = "bold"))`-42.51993430*exp(-10002901.82*t)+0.1702286306e-1*exp(-0.1544864054e11*t)*`#msub(mi("V"),mo("in",fontweight = "bold"))`+.43114634*exp(-0.1544864054e11*t)+24., v[gs](t) = -0.3674479877e-3*exp((30764029.75+1577330139.*I)*t)*`#msub(mi("V"),mo("in",fontweight = "bold"))`-0.8037470713e-2*exp((30764029.75+1577330139.*I)*t)-(0.7193539082e-1*I)*exp((30764029.75+1577330139.*I)*t)-(0.3899911367e-2*I)*exp((30764029.75+1577330139.*I)*t)*`#msub(mi("V"),mo("in",fontweight = "bold"))`-0.6005358003e-4*exp(-10002901.82*t)*`#msub(mi("V"),mo("in",fontweight = "bold"))`-8.504088823*exp(-10002901.82*t)+0.1556340602e-2*exp(-0.1544864054e11*t)*`#msub(mi("V"),mo("in",fontweight = "bold"))`+0.3941819611e-1*exp(-0.1544864054e11*t)+12.}

The <maple> tag doesn't recognize the typesetting codes for V_in: make that

{v[ds](t) = .4958026528*exp((30764029.75+1577330139.*I)*t)*V[`in`]+9.149694069*exp((30764029.75+1577330139.*I)*t)-(.766755745*I)*exp((30764029.75+1577330139.*I)*t)-(0.3290893603e-1*I)*exp((30764029.75+1577330139.*I)*t)*V[`in`]-0.3002643e-3*exp(-10002901.82*t)*V[`in`]-42.51993430*exp(-10002901.82*t)+0.1702286306e-1*exp(-0.1544864054e11*t)*V[`in`]+.43114634*exp(-0.1544864054e11*t)+24., v[gs](t) = -0.3674479877e-3*exp((30764029.75+1577330139.*I)*t)*V[`in`]-0.8037470713e-2*exp((30764029.75+1577330139.*I)*t)-(0.7193539082e-1*I)*exp((30764029.75+1577330139.*I)*t)-(0.3899911367e-2*I)*exp((30764029.75+1577330139.*I)*t)*V[`in`]-0.6005358003e-4*exp(-10002901.82*t)*V[`in`]-8.504088823*exp(-10002901.82*t)+0.1556340602e-2*exp(-0.1544864054e11*t)*V[`in`]+0.3941819611e-1*exp(-0.1544864054e11*t)+12.}

You could also use evalc to express the complex exponentials in terms of real exponentials times sines and cosines.

Note that squaring a complex quantity doesn't make it real.  Perhaps you mean to square the absolute value.  Again, you probably want to use implicitplot3d to plot a surface abs(u(r,theta,z)) = c for some suitable value of the constant c.

Note that squaring a complex quantity doesn't make it real.  Perhaps you mean to square the absolute value.  Again, you probably want to use implicitplot3d to plot a surface abs(u(r,theta,z)) = c for some suitable value of the constant c.

You might try using Classic rather than Standard GUI, and exporting programmatically using plotsetup rather than from the menu.  Thus:

> P:= plot3d(whatever):
  plotsetup(eps, plotoutput="c:/mypath/myplot.eps",plotoptions=
   "portrait,noborder,height=3in,width=3in,leftmargin=0in,bottommargin=0in");
  P;
  plotsetup(default);

You might try using Classic rather than Standard GUI, and exporting programmatically using plotsetup rather than from the menu.  Thus:

> P:= plot3d(whatever):
  plotsetup(eps, plotoutput="c:/mypath/myplot.eps",plotoptions=
   "portrait,noborder,height=3in,width=3in,leftmargin=0in,bottommargin=0in");
  P;
  plotsetup(default);

First of all, make sure the function you're trying to plot will give numerical values when given numbers for r, theta, z.  If you gave us the actual Maple input, we might be able to help with that.

Second, what sort of a plot are you looking for?  You have a function of three variables.  It should give you a value at each point of three-dimensional space. 
It's impossible to plot that directly (the graph would have to be a four-dimensional object).  You might, for example, plot the surface f(r, theta, z) = c for some constant c.  This could be done using implicitplot3d in the plots package.  For example:

> plots[implicitplot3d](r^2*sin(theta)^2*exp(-r^2-z^2) = 1/10,
   r = 0 .. 2, theta = -Pi .. Pi, z = -2 .. 2,  
   coords=cylindrical,  view=[-2..2,-2..2,-2..2],  
   axes=box, labels=[x,y,z], style=patchcontour,
   grid=[40,40,40],orientation=[30,60]);

Note: because of a bug, this doesn't work properly without the view option.

First of all, make sure the function you're trying to plot will give numerical values when given numbers for r, theta, z.  If you gave us the actual Maple input, we might be able to help with that.

Second, what sort of a plot are you looking for?  You have a function of three variables.  It should give you a value at each point of three-dimensional space. 
It's impossible to plot that directly (the graph would have to be a four-dimensional object).  You might, for example, plot the surface f(r, theta, z) = c for some constant c.  This could be done using implicitplot3d in the plots package.  For example:

> plots[implicitplot3d](r^2*sin(theta)^2*exp(-r^2-z^2) = 1/10,
   r = 0 .. 2, theta = -Pi .. Pi, z = -2 .. 2,  
   coords=cylindrical,  view=[-2..2,-2..2,-2..2],  
   axes=box, labels=[x,y,z], style=patchcontour,
   grid=[40,40,40],orientation=[30,60]);

Note: because of a bug, this doesn't work properly without the view option.

It works perfectly for me: I get

  0., 0., 0., 0., 0., 0.

which is correct: cos(0)*cosh(0) - 1 = 0.

Perhaps you want different solutions.  You could specify intervals for beta.  For
example:

> seq(fsolve(cos(beta*L)*cosh(beta*L)-1=0, beta = Pi/L .. 2*Pi/L),
   L = 1 .. 6);
 

4.730040745, 2.365020372, 1.576680248, 1.182510186, .9460081490, .7883401241

It works perfectly for me: I get

  0., 0., 0., 0., 0., 0.

which is correct: cos(0)*cosh(0) - 1 = 0.

Perhaps you want different solutions.  You could specify intervals for beta.  For
example:

> seq(fsolve(cos(beta*L)*cosh(beta*L)-1=0, beta = Pi/L .. 2*Pi/L),
   L = 1 .. 6);
 

4.730040745, 2.365020372, 1.576680248, 1.182510186, .9460081490, .7883401241

1) Check your spelling. 

2) Now you have n, m and q as input parameters, but you're not using them as such. 
Instead you're using them as loop indices.  Never do that. Loop indices should be local variables.

3) Your finalmatrix[n,m]; does nothing.

4) You need to assign an entry of finalmatrix after finishing the loop on q, not inside it.

5) You should check whether c1 = r2 before starting your loops, not within them. 
c1 <> r2 should result in an error. not just printing "not calculatable".  You can use an error statement to do this.

6)  Your return statement makes no sense, and shouldn't be inside the loops.  You want to return finalmatrix at the end of the procedure, not in the middle.

1) Check your spelling. 

2) Now you have n, m and q as input parameters, but you're not using them as such. 
Instead you're using them as loop indices.  Never do that. Loop indices should be local variables.

3) Your finalmatrix[n,m]; does nothing.

4) You need to assign an entry of finalmatrix after finishing the loop on q, not inside it.

5) You should check whether c1 = r2 before starting your loops, not within them. 
c1 <> r2 should result in an error. not just printing "not calculatable".  You can use an error statement to do this.

6)  Your return statement makes no sense, and shouldn't be inside the loops.  You want to return finalmatrix at the end of the procedure, not in the middle.

Like this:
 

> subs(RGB=HSV,plot3d(0, x=-3..3,y=-3..3,colour=[argument(-ln(x+I*y)),-abs(ln(x+I*y)),1],
    axes=box,orientation=[-90,0],grid=[101,101],style=patchnogrid));

Like this:
 

> subs(RGB=HSV,plot3d(0, x=-3..3,y=-3..3,colour=[argument(-ln(x+I*y)),-abs(ln(x+I*y)),1],
    axes=box,orientation=[-90,0],grid=[101,101],style=patchnogrid));

If I1, I2, I3 and I4 are "knowns", you shouldn't assign values to them as you did.  I think you want something like this.  I'll call the values i1, i2, i3, i4.

> solve({I1=i1,I2=i2,I3=i3,I4=i4},{M,a,b,theta},AllSolutions);

Unfortunately this doesn't produce a result, so I'll try Groebner basis methods.

> with(Groebner);
G:= map(numer,subs(cos(theta)=c,sin(theta)=s,convert([I1-i1,I2-i2,I3-i3,I4-i4,s^2+c^2-1],trig)));
B:= Basis(G,plex(a,b,c,s,M,i1,i2,i3,i4));
B[1];

i4*i3^2*i2^2-i3^2*i2^2*i1-i2^2*i1*M^2+i2*i1^2*M^2-i4^2*i3*i1^2+i4^2*i2*i1^2+i4*i3^2*M^2+i3*i2^2*M^2*I+i4^2*i1*M^2*I+i4^2*i3^2*i1*I-I*i4*i2^2*i1^2-I*i4*i1^2*M^2-I*i4^2*i3^2*i2-I*i3^2*i2*M^2+i3*i2^2*i1^2*I-i4^2*i3*M^2+(1-I)*i4*i3^2*i2*i1+(-1+I)*i4*i3*i2^2*i1+(1-I)*i4*i3*i2*i1^2+(-2+2*I)*i4*i3^2*i2*M-(1+I)*i4*i3*i2^2*M-(1+I)*i4*i3^2*i1*M+(1+I)*i3^2*i2*i1*M+(1+I)*i4*i2^2*i1*M+(2-2*I)*i3*i2^2*i1*M+(1+I)*i4*i3*i1^2*M+(-2+2*I)*i4*i2*i1^2*M-(1+I)*i3*i2*i1^2*M+(1-I)*i4*i3*i2*M^2+(-1+I)*i4*i3*i1*M^2+(1-I)*i4*i2*i1*M^2+(-1+I)*i3*i2*i1*M^2+(-1+I)*i4^2*i3*i2*i1+(1+I)*i4^2*i3*i2*M+(2-2*I)*i4^2*i3*i1*M-(1+I)*i4^2*i2*i1*M

This isquadratic in M, and has two solutions in term of i1,i2,i3,i4 (which I won't show here because they're rather complicated).

> solve(%,M);