that you clean up your supplied worksheet. I started to go through it, removing (apparenlty?) redundant code, erroneous code etc in the hope that I might actually understand what you are trying to achieve

But there is just so much irrepevant stuff in your worksheet, that frankly, I got a bit lost. I suggest that you read the attached very carefully and construct a worksheet with a minimum of extraneous rubbish

badCode.mw

The attached shows how to use dsolve() with the output=Array() option. This will provide results for all unknowns at the values specified in the Array definition. As a test I have overlaid plots from dsolve() with and without the output=Array() option - and they coincide exactly

odeProb.mw

You state that

Yours, Tomleslie, works fine. Instead of getting some 5 columns indexed (t,x(t),vx(t),y(t),vy(t)) I get some 7 columns one containing ax(t) and ay(t). But if I don't do exactly what you write, it does'nt work anymore. For instance if I define the dummies before the eqs, or if, instead of defining the dummies, I onely write something like :

Ax:=t->diff(x(t),t,t);Ay:=t->diff(y(t),t,t);

And I don't understande why it does'nt work anymore.

One more time, thank you both.

Unless you actually upload the code which is non-functional, I have no idea how to diagnose your problem. Use the big green up-arrow in the MaplePrimes toolbar to upload code which is doing something you don't want/expect and I will do my best to help

But sometimes alternative approaches are good

Basically, the dsolve(,,numeric) process does not return explicit expressions for the highest order derivatives in the supplied ODE system. You can argue that maybe it should -but if you think about it, if all other lower-order derivatives and function values are computed, then you can compute the highest-order derivatives by substitution in the original equation(s). This is esentailly what Rouben has done - and it works!! I have no argument with this at all.

An alternative approach (I think(?) first posted here by Preben Alsthom) is to introduce "dummy" variables, setting these equal to the highest derivatives in the system. This would mean rewriting Rouben's example as

restart;
de1 := diff(x(t),t,t) = -x(t) - diff(x(t),t) + y(t);
de2 := diff(y(t),t,t) = -y(t) - diff(y(t),t) + x(t);
dummy1:=Z1(t)=diff(x(t),t,t); # introduce dummy variable Z1(t)
dummy2:=Z2(t)=diff(y(t),t,t); # introduce dummy variable Z2(t)
ic := x(0)=0, D(x)(0)=1, y(0)=0, D(y)(0)=1;
dsol := dsolve({de1, de2, ic, dummy1,dummy2 }, numeric);

and then, when you want to produce plots of the highest derivatives, all you need is

plots:-odeplot(dsol,[t,Z1(t)], t=0..10);
plots:-odeplot(dsol,[t,Z2(t)], t=0..10);

we are still living with the consequences of an era when those programming GUIs chose to use fixed pixel size, rather percentage of screen size as a display unit. As pixel sizes got smaller (aka UHD monitors), then GUI entities became progressively less readable.

You may (or may not!) be surprised just how many GUIs of how many software packages are affected by this

Only way I know to get around this, is rto educe your display resolution - easy on Windows, just right-click on the desktop, select Display Resolution - and keep reducing until everything becomes more "visible", if a little "fuzzy".

I assume that Linux/Mint has a similar capability to reduce screen resolution

If, as you say, the solution is given as

x_k = x_k for all variables

then your system is solved for any values of your variables. So there are an infinite number of solutions. I'm guessing that this isn't what you want

Consider the procedure

f := proc(x, rec := true)
          if   rec
          then Optimization:-NLPSolve( 2,
                                       y -> (y -~ 1/4).(y -~ 1/4),
                                       [],
                                       [<-1, -1>, <1, 1>],
                                       maximize = false
                                     )
          end if;
          (x -~ 1/4).(x -~ 1/4)
     end proc:

This will only ever return the result of the final 'command', ie (x -~ 1/4).(x -~ 1/4). The Optimization() command within this procedure may or may not be executed, depending on

1. if the procedure is called without a second argument, the defualt rec:=true means then the Optimization() command will be executed.
2. if the procedure is called with a second argument equal to 'true' then the Optimization() will be executed
3. if the procedure is called with a second argument equal to 'false' the the Optimization() command will not be executed

However it really does not matter whether the Optimization() command is executed or not - since its output is never used for anything and will never be returned by this procedure - the result of the final command, ie (x -~ 1/4).(x -~ 1/4) will always be returned!

Irrespective of input arguments, and since the outcome of the Optimization() is never used/returned/orAnything, then the above procedure can be written much more efficiently as

f := proc(x, rec := true)
              (x -~ 1/4).(x -~ 1/4)
      end proc:

and you don't even need a default for the second argument - becuase obviously it is also irrelevant

You appear to have some kind of weird typographic problem in one of the terms of the set 'bcs3'. It *looks* like v/(1-v^2), and *should* evaluate accordingly - but it doesn't! Maple is treating it as a 'name'. Since this 'name' does not evaluate to anything, Maple considers it as an "unknown". The presence of this extra "unknown" is why Maple "thinks" it needs an extra boundary condition in the dsolve() process.

I managed to isolate/retype the offending term in bcs3. (Trust me, since conventional find/replace does not work well, if you choose to use 2D input, this is a lot harder than it looks).

The attached worksheet now executes. However as you loop through the "extra_bcs", some of the dsolve() commands give accuracy warnings, and the occasional one fails with "Matrix is singular" errors. I have made no attempt to address any of these later issues

pdeProb.mw

Your code could do with a serious "tidy-up"

The "limiting number of evaluations reached" warning can be overcome by setting the evaluationlimit option in nteh Maximize command. I set this to 1000 in the attahced and both of your examples run in negligible time.

You claim there is a "gap" in one of these plots - I agree there is a gap. Can only think that this occurs as a result of you algorithm, in particaular your use of piecewise functions, where continiuity will never be "guaranteed"

Check the attached

Download lev.mw

Have to confess that I'm surprised that this

restart;
expr:=4/x^(1/3)+1=0;
solve(expr, x);

does not work, particularly when this

restart;
expr:=4/x^(1/3)+1=0;
isolate(expr, x);

provides the answer you want.

The only(?) way I can persuade solve() to return the required value is to use RealDomain(), as in

restart;
with(RealDomain):
expr:=4/x^(1/3)+1=0;
solve(expr);

you can use simplify(..,power), or by appropriate use of freeze()..thaw(), you can simplify only some powers, as in

   restart;
   expr:= a^k*k/(a*a__0^k);
#
# either ("simplifies" all powers)
#
   simplify(expr,power);
#
# or (don't simplify powers of a__0)
#
   thaw
   ( simplify
     ( subs
       ( a__0^k=freeze(a__0^k),
         expr
       ),
       power
     )
   );

Using only the commands with the geom3d package, use of the solve() command is unnecessary. Your objective can be achieved with the code (using only commands from the geom3D package)

   restart;
   with(geom3d):
#
# Define the plane
#
   plane( p, 2*x+y-3*z-3 = 0, [x, y, z]):
#
# Define the line
#
   line( l, [ -1+3*t,1+t, 2-t ],t ):
#
# Draw the line and plane, showing the
# intersection
#
   draw([p, l]);
#
# Compute the intersection
#
   intersection(Z, l, p):
#
# Get the details of the intersection point
#
   detail(Z);

#
# If you only want the coordinates of the intersection
# point, then the following will work
#
    coordinates(Z);

One warning - don't call the point of intersection between the plane and the line 'I'. This gives me a slightly odd error message . I suspect that the problem is that (by default), within Maple 'I' refers to square root of -1.

If you want the short version, without any explanations/comments, then the following will also supply txhe answer you require

 restart;
 with(geom3d):
 plane( p, 2*x+y-3*z-3 = 0, [x, y, z]):
 line( l, [ -1+3*t,1+t, 2-t ],t ):
 coordinates(intersection(Z, l, p));

two of which are too long to present here, but the third is *close* to that , presented by the OP

Download sysSol.mw

Just because alternatives are good,

restart;
with(plots):
f:= w^8+w^6+4*w^4+w^2+1:
g:= w^16+2*w^14+9*w^12-2*w^10+44*w^8-2*w^6+9*w^4+2*w^2+1:
p1:=complexplot([fsolve(f, complex)], style=point, symbol=solidcircle, symbolsize=20, color=red):
p2:=complexplot([fsolve(g, complex)], style=point, symbol=solidcircle, symbolsize=20, color=blue):
p3:=complexplot(cos(theta)+I*sin(theta), theta=-Pi..Pi, scaling=constrained, color=green):
display([p1,p2,p3]);

In Maple 2016, if I access the help page for Feynman diagrams using ?FeynmanDiagrams. then, within  the help window use View->Open Page as Worksheet, I obtain a worksheet, which executes (using !!!!) with no errors/problems

In Maple 2017, if I repeat the same process,the resultant worksheet contins most  of the errors you llist - in other words something got broken.

This is worthy of an SCR

As has been observed, it is true that Maple's Tolerances package only handles symmetric tolerances.

However, reading the help pages for the Tolerances package, indicates that this is more-or-less a "convenient interface" to the evalr() command, which evaluates expressions over "ranges" - known, for the purposes of that command as INTERVAL(s). Rather than using the Tolerances package, it is not too difficult to set up something which uses the evalr() command directly, and which will permit the use of asymmetric ranges.

I (personally) find that the amount of typing required to use evalr() stuff is excessive, so the attached contains some "helper" function definitions, which allow me/you to minimise typing, and handle arbitrary amounts of variables and ranges.

One drawback/annoyance is how to conveniently express a parameter with asymmetric tolerances. I have chosen to do this as a list of three values given by parameter=[nominalValue, negTol, posTol], so that actual values for the parameter occupy the range nominalValue+negTol..nominalValue+posTol. The solution is returned in the same format

Can you "break" the attached with some combination of parameterValues, parameterTolerances, functionDefinitions - almost certainly. Otherwise I would have managed to produce something (in a few lines) which is even better than Maple's Tolerances package - which is exceedingly unlikely.

So with all of thse caveats, check out the attached

Download tolProb.mw