I tried in both Maple 16 & 17, obtaining the same results. The "hang" occurs in simplify. I don't know whether it is an infinite loop, or just an expression that is too complicated for simplify. Since the expression fits on one line on my screen, the latter seems unlikely. Running with infolevel[simplify]:= 5, the only information is a seemingly endless repetition of

• simplify/do: applying  commonpow  function to expression
• simplify/do: applying   power  function to expression

The hang will occur for any of z[17], z[18], z[19], z[20]. It does not matter whether the expand or factor or both are applied before the simplify. Only the first two solutions need to be substituted in for the hang to occur. So here is one of the seemingly simple expressions for which it hangs:

subs(solution1, solution2, z[20]);

lprint(%);

0 = -(-b4*f8+b8*f4)*conjugate(f7)*d5*(conjugate(b4)*conjugate(d8)-conjugate(b8)*conjugate(d4))/((a4*b8-a8*b4)*(conjugate(a4*b8)-conjugate(a8*b4)))-(-a4*f8+a8*f4)*conjugate(f7)*d5*(conjugate(a4)*conjugate(d8)-conjugate(a8)*conjugate(d4))/((a4*b8-a8*b4)*(conjugate(a4*b8)-conjugate(a8*b4)))+conjugate(e7)*e5

How about diff(z(t), t) = m with initial condition z(0) = n? Then the system becomes

diff(x(t), t) = -x(t) + 2*z(t)*y(t) - 2*k*z(t)^2,

diff(y(t), t) = -y(t) + k*m - z(t)*(x(t)/2 - 1),

diff(z(t), t) = m;

Here it is, by guiding Maple at each step (or "holding its hand").

Download int_sec^n.mw

You need to create an Array (or Matrix) before you can assign to it. Otherwise, the object becomes a table.

(n,m):= (3,3);

U:= Array(0..n+2, 0..m):

U[.., 0]:= < y(a), seq(f[k], k= 0..n), y(b) >;

Note that it is U[.., 0] rather than U[0, ..]. The latter would create the first row rather than the first column.

First, you must create and store S as a Vector, list, or table---not a set. There are two problems with using a set:

1. Duplicate entries are eliminated, and in your case there is a small chance that two of the numbers associated with an index pair are the same.
2. A set is stored sorted, and you'll want to avoid the super-linear overhead of the sort.

Which one is the most efficient to use---Vector, list, or table---depends entirely on the code used to create S. So, I'd like to see that code. Note that if you use a table, then the sum can be accumulated as the entries are created. Create it with table(sparse) so that entries are initialized to 0.

The process of converting S to a Matrix is independent of which of the three formats that you use. I will assume that the Matrix order N is known a priori; if instead the order needs to be determined from the scan of S, then my code will need to be modified.

Download sparse.mw

Because of the issue of the target lists possibly appearing as sublists, I recommend that you use strings instead of lists as your data structure. Then you can use StringTools:-SubstituteAll for the operation.

Even if you need to use lists, it may be faster to perform the operation by converting to a string and then back to a list every time. That's because the StringTools commands are written in compiled code.

First, get rid of linalg. If you need linear algebra (you don't for this bit of code), use LinearAlgebra instead.

Only the terminal punctuation of the whole for-do loop, either colon or semicolon, makes any difference in the printed output. So we need to suppress all the output by ending the loop with a colon, and, within the loop, selectively printing what we want with a print statement.

for h from .5 by .1 to 1 do
     eq1 := x*(-b*x^2-x+1);
     eq2 := y*((a*x*x)/(b*y^2)-d-h*y);
     S := solve({eq1, eq2}, {x, y});
     SS := solve(subs(S[3], {omega^4+(h*y+x)*omega^2+h^2*x-y}), {omega});
     tau := simplify(subs(S[3], subs(SS[3], (b^2*h*y+a*x)/omega)));
     print('tau' = tau);
end do:

You wrote:

what if we change the original equation to a pde with the derivative on the left hand side to a derivative with repect to t?

Then the rest of your code in your original question works. The animation shows a very interesting wave pattern. Thus I think this is what the equation was origianlly intended to be. It doesn't make sense to me to do it in 3D though.

S:= pdsolve(pde, {bc, ic}, numeric):
S:-animate(t= 0.1, frames= 64);

Your matrix is severely rank deficient, a 57x57 with rank 44. The matrix is singular. I don't think it's a problem in your program; rather, it's a problem with your initial data. You could seek a least-squares solution to the linear system. However, I know nothing about this type of engineering problem, so I don't know if that's appropriate.

You should not use the package linalg. Use LinearAlgebra instead.

By the way, this is the most primitive Maple code of its size that I have ever seen---although it is neatly presented. It's very easy to make a mistake if you initialize huge matrices by indexing each element in its own assignment statement.

You wrote:

can this function be animated in 3d?

Sure. It's very easy:

S:= pdsolve({pde, bc}):
plots:-animate(plot3d, [eval(U, S), x= 0..1, t= 0..T], T= 0..10);

This particular function makes a rather boring animation though.

Why are you using Vectors instead of lists? You can use Vectors, but for every operation of this type, you'll need to convert to a list and then back to a Vector. And possibly you'll need to convert to and from a set also.

Once you have it in list form, there are several commands in ?ListTools to help with the things that you're asking about, specifically, ?ListTools,Categorize ?ListTools,Classify ?ListTools,FindRepetitions ?ListTools,MakeUnique and ?ListTools,SearchAll

You do it by assigning the output of fsolve to a variable (let's say Sol), and then using Sol as the second argument to eval with the first argument being the expression which contains the variables that fsolve solved for. (This same technique can be applied to any of the "solve" commands: dsolve, solve, rsolve, pdsolve, isolve, msolve, etc.)

Example:

Sys:= {'randpoly([x||(1..3)])+randpoly(x, degree= 0)' $ 3}; 
 /           3        3   2       3           3              3
{ 31 x1 x2 x3  - 27 x2  x3  - 2 x1  x2 - 31 x2  x3 + 25 x2 x3 
 \                                                            

                            4        5       4                 2
   - 97 x2 x3 + 65, 68 x1 x2  - 29 x2  + 5 x2  x3 + 73 x1 x2 x3 

             3           2            5        4   
   + 95 x1 x3  + 31 x1 x3  - 26, 91 x1  - 22 x1  x2

                3        2   2        4                   \ 
   + 51 x1 x2 x3  + 52 x1  x3  + 36 x3  + 18 x1 x2 x3 - 27 }
                                                          / 
Sol:= fsolve(Sys, {x||(1..3)});
      {x1 = -0.194514327765446, x2 = -0.806303974671123, 
        x3 = -0.998794548580666}

plot(eval(x1*t^2 + x2*t + x3, Sol), t= -2..2);

To add static frames to an animation, you first create the animation (using, for example, plots:-animate or plots:-display(..., insequence= true)), and then you use display again to merge the animation with a single static plot, which becomes a fixed part of every frame. This merged animation can be played as is, or it can be merged with another animation, once again with didplay. There's no limit to how many levels of merging with display that you can do. So, you can create animations that have parts that are static during only part of the animation.

So, here is your code with my modifications. It uses four levels of display, whereas your original had two levels:

restart: 
with(plots): 
auto:= (x, y, C)-> plots:-pointplot([[y, x]], color= C):

p1:= animate(auto, [t, .17+0.55e-1*t, blue], t= .2 .. -2, frames= 10):
p2:= animate(auto, [t, .16+0.45e-1*t, red], t= .7 .. -2, frames= 20):
p3:= animate(auto, [t, .15+0.45e-1*t, green], t= 1.2 .. -2, frames= 30):
p4:= animate(auto, [t, .15+0.45e-1*t, black], t= 1.7 .. -2, frames= 40):

p5:= animate(auto, [t, .12+0.45e-1*t, blue], t= 2.2 .. .2, frames= 50):
lastframe5:= auto(eval([t, .12+0.45e-1*t, blue], t= .2)[]):

stoppt:= frame-> 2.7+(.7-2.7)*frame/60: #Split t-range for auto 6.
p61:= animate(auto, [t, .11+0.45e-1*t, red], t= 2.7 .. stoppt(50), frames= 50):
#Frames 51-60 for auto 6:
p62:= animate(auto, [t, .11+0.45e-1*t, red], t= stoppt(51) .. .7, frames= 10):
#Merge animations for auto 6:
p6:= display([p61, display([p62, lastframe5])], insequence):
lastframe6:= auto(eval([t, .11+0.45e-1*t, red], t= .7)[]):

firstframe7:= auto(eval([-.17+0.35e-1*t, t, green], t= -.2)[]):
p7:= animate(auto, [-.17+0.35e-1*t, t, green], t= -.2 .. 7, frames= 20):

firstframe8:= auto(eval([-.15+0.25e-1*t, t, black], t= -.7)[]):
p8:= animate(auto, [-.15+0.25e-1*t, t, black], t= -.7 .. 6, frames= 20):

A:= display([firstframe||(7..8), p||(1..6)]):
B:= display([p||(7..8), lastframe||(5..6)]):
display(
   [A, B], insequence, 
   scaling= constrained, symbol= solidbox, symbolsize= 20,
   title= ``
);

The improved animation:

Your original animation:

I have a feeling that this can be handled by the events option to dsolve (see ?dsolve,events ), but I don't know the details. There are far too few examples on the help page and an overwhelming number of options described above the examples.

But here is a quick-and-dirty (linear search) procedure that will give you the last values plotted by odeplot (the last X value and the last Y value) before the first undefined values. (I'd be pleased if someone else could comment on a better way to do this.)

LastValues:= proc(P::specfunc(anything,PLOT))
     local M,X,Y,m,n;
     M:= op([1,1], P);
     X:= convert(M[.., 1], list);
     m:= ListTools:-Search(HFloat(undefined), X);
     Y:= convert(M[.., 2], list);
     n:= ListTools:-Search(HFloat(undefined), Y);
     `if`(m < 2, undefined, X[m-1]), `if`(n < 2, undefined, Y[n-1])
end proc;

Call the procedure as, for example, LastValues(plott[1]), and it should return two numbers, the first being the last X value and the second being the last Y value.

Please let me know if this does what you need.

dotprod(j,k) is 0 after you radnormal or simplify it. There is nothing wrong with your code.

You should not use linalg anymore. Use LinearAlgbera. It is available in Maple 8.