When defining the s-variable for transfer functions s=I*omega, so the following

restart;
er:=(0.1299e-4*(17.1740*Omega^2-1.570200000*10^6))*Omega^2/(-196.1270800*Omega^4+3.954121290*10^7*Omega^2-1.877174100*10^12);
er := subs(Omega=s/I, er):
with(DynamicSystems):
sys := TransferFunction(er):
BodePlot(sys, range=200..500, linearfreq=true, numpoints=1000);

get pretty close to the textbook plot. You might be able to get closer by fiddling witht the plot option, scaling etc

Consider your second "obvious" answer, x=-I. If I substitute this in your expression, I get abs(I+ -I) which is abs(0), which you appear to think is equal to 1 - on which planet?

Consider now the two solutions which are returned x=1-I and x=-i-I, Well

abs(I+x)=abs(I+(1-I))=abs(1)=1 so this is correct
abs(I+x)=abs(I+(-1-I))=abs(-1)=1, so this is correct

The "interesting" solution which you propose is x=0, so that

abs(I+x)=abs(I+0)=abs(I), and since this is equal to 1, I'm not sure why solve() didn't find it :-(

The three "correct" solutions would eem to be -1-I, I, i-I

For 2D plots, one of the easiest ways is to use the plot option 'colorscheme'. How 'pretty' you want to make this look is up to you, but consider the simple pointplot, with the simple 'linear' colorscheme given by

restart;
d1:=[seq( [x, 1/x], x =1..100)]:
plots[pointplot](d1, colorscheme=[red, yellow, blue]);

And there are many other ways to define the 'colorscheme' option - see ?colorscheme and experiment till you get what you want!

A better bet than Riemann sum might be to use the evalf(Int()) approach. See the help page at ?Numerical Integration.

This allows the user to set a variety of methods, including several with 'quadrature' in their names. Since you do not supply the actual integral, I can only suggest that you experiment with these

Somehow(?) you have managed to turn on numeric formatting for one or more execution groups. This option is normally accessed from the toolbar menu Format->Numeric Formatting.

In theory if you turn on some kind of numeric formatting then it applies to the current execution group and all subsequent groups, so the latter ought(??) to handle the format.

On the other hand it is possible to turn it on for (say) one execution group, then off for the next - so the formatting only applies for a single group. Then copying from the group with numeric formatting set, to anywhere is probably going to cause problems

I notice that you are using some aspects of the Units package. I don't *think* that any command/option within this package should change the numeric formatting anywhere but I'm not 100% sure......

If you simply ask Maple to solve this equation (without specifying RealDomain), then

solve(sqrt(x)+sqrt(-x^2+1) = sqrt(-4*x^2-3*x+2), x)

will return

-1-2^(1/2), -5/9+(1/9)*34^(1/2)

ie, exactly the same as Mathematica.

When you invoke the RealDomain package in Maple, then as the help page tells you

By default, Maple performs computations under the assumption that the underlying number system is the complex field. The RealDomain package provides an environment in which computations are performed under the assumption that the basic underlying number system is the field of real numbers.

Now given that Maple is performing computation on the basis that the underlying number system is the field of real numbers, what would you expect Maple to do with the first term of your  expression and the first "solution", ie sqrt(-1-2^(1/2)). This quantity is not real.

Just to complicate things slightly, when maple 'parses' an input expression, then it does some kind of 'basic simplification', and sometimes this can make complex quantitie " disappear". For example

RealDomain:-solve( sqrt(-1)^3+x=1-I,x)

will return 1 (correctly), becuase the initial simplification process will eliminate the complex quantities, before being passed to the RealDomain[solve]()

Well a certain amount depends on whether you want to "zoom" all three dimensions by the same amount!

The attached shows examples for a couple of 'toy' functions, where

1. all (ie x,y,z) directions are "zoomed" by the same amount,and
2. where x,y dimesions are defined and only the z-direction is "zoomed"

plotZoom.mw

Replacing the "toy" functions with anything you want should be easy, although I suggest avoiding anything which goes infinite over the range you desire, because this will screw the calculation of the range in 'z', and hence the 'zoom' in 'z'.

I could deal with this case if I had to, but if I don't have to.......

1. You could try reading the manual, or
2. You could try Vector( 101, i->(i-1)*10);

Well if I do a quick check using

  restart;
#
# Get Maple solution from dsolve()
#
  eqn:= diff(y(t),t$3)-12*diff(y(t),t$2)+48*diff(y(t),t)-64*y(t)=12-32*exp(-8*t)+2*exp(4*t);
  mapleSol:=dsolve(eqn);
#
# Type in "hand" solution *very carefully*
# NB I have changed c1, c2, c3 in OP's
# original to _C1, _C2, _c3 respectively
#
  handSol:=y(t)=_C1*exp(4*t)+_C2*t*exp(4*t)+_C3*t^2*exp(4*t)-3/16+(1/54)*exp(-8*t)+(1/3)*t^3*exp(4*t);
#
# Check whether handSol=mapleSol
#
  is(simplify(mapleSol)=simplify(handSol));

The output of the final command is 'true', indicating that the Maple solution is equal to the "hand" solution

PS This probably ought to be in the MaplePrimes/Questions section, rather than MaplePrimes/Posts

1. You want a recursive definition of p(t) in terms of p(t-1).
2. In order to achieve this there has to be some kind of "start" condition, otherwise the recursion is infinite.
3. I have (arbitrarily) decided  that p(-1)=0, so that only arguments >=0 can be used.

Given the above, the following will work

restart;
#
# Create recursive function definition
#
   p:=x-> `if`( x=0,
                      m[tau]/(1-psi*(beta-1)),
                      m[tau]/(1-psi*(beta-1)) - psi*(beta-1)*p(x-1)/(1-psi*(beta-1))
                    );
#
# compute/assign the first 6 entries for the
# table 't'
#
   seq( assign(t[j], p(j) ), j=0..5);

There are more "efficient" ways to do this (using a recursive procedure with a 'remember' table rather the simple function above), but unless you are going to calculate an *lot* of terms, you probably won't notice the difference

I have not attempted to make any simplifications oon the expressions which are generated as table entries

Seems like you have a very complicated way to achieve what seem to want.

You generate tables, convert to lists, then vectors. This isn't wrong but it does make things unnecessily complicated.

So far as I can tell, the following generates, fits, and plots, the same data

restart;
with(Statistics):
loopStart:=0:
loopEnd:=2:
loopStep:=0.01:
sz:=loopEnd/loopStep:
myt:= Array( 0..sz, i->loopStep*i):
myA:= Array( 0..sz, fill=0.):
for j from 1 by 1 to sz do
    myA[j]:=myA[j-1]+0.01*(loopStep*j)^2;
od:
p1:=plot( myt, myA, style=point, color = red):
poly1:=PolynomialFit(10, myt, myA, time):
p2:=plot(poly1, time=0..2, color=blue):
plots[display]([p1,p2]);

NB I don't think you need a 10-th order polynomial to fit a curve this smooth.

Your outer loop contains the statement

q:=rhs(sol[3]);

The second time you execute the outer loop, the statement

ics:=q(0)=1, p(0)=n;

will therefore be equivalent to

ics:=rhs(sol[3])(0)=1, p(0)=n;

and since sol[3] was also assigned in the first loop iteration, then rhs(sol[3])(0) will evaluate, and return 1.0, so that this is equivalent to

ics:=1.0=1, p(0)=n;

which is of course meaningless as boundary condition.

'Minimal change' fix is probably to add the statement

unassign(q);

immediately prior to the 'end do' statement of the outer loop. Loop will now execute - does for me anyway.

NB Suggest you change the semicolons on the 'end do' statements to colons - otherwise you are going to generate an awful lot of verbiage you don't really need

Answers I obtained for  F, F1, S, S1, T, T1 are 247, 500, 271, 500, 1322, 4000 respectively. Obviously I have no idea whether these are correct (or what you expect), since I have not reallly anysed what you are trying to achieve

So far as I can tell, fsolve() is unable to find a solution, so the command returns unevaluated, which breaks all the subsequent code.

The DirectSearch package did return a solution, namely

A = -.6839299669509255427872547,
B = -.7058228818317920507398118,
C = .2983804370023421519862147,
E = -1.000000000023554794609445

Residuals are a little on the high side, though

If I ignore all the superfluous statement, you problem is with the statement

bcs2 := eval[recurse](convert(ds[1], D), `union`({y = 0}, bcs3));

In particular with the part y=0. Before or after the conversion of ds[1], it contains several terms whose denominator is either y, y², y³,or even y⁴.

Now with y=0, what do you expect to happen to these terms in ds[1]?

N:=infinity;

ought to work.

I have just tried it - know what, it does work!