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.

I can't tell if you want to create a surface, or if you simply want the sequence of curves in a 3D plot. If you just want the curves, it requires merely a trivial modification to your code:

1. Change the curve that you plot in the odeplot command from [X(t), Y(t)] to [ j, Y(t), X(t)].
2. Change the display command to display(convert(plott, list)).
3. (Optional) Include an option such as thickness= 5 in the odeplot to improve visibility.

If you do want to create a surface, let me know.

[I converted your Post into a Question.-- Carl Love as a moderator]

Please upload a worksheet and/or post your equations in plaintext so that we do not need to retype them into Maple. Thanks.  The solution to your problem is probably to simply include a range in the fsolve command. See ?fsolve,details

fsolve({f, g}, {x= 0..Pi/2, y= 0..Pi/2});

Here's how to do it with a continuous transformation to your existing color function, which is presumed to return a value between 0 and 1 (the HUE color scale). Keeping it continuous is very very nice when you want colors to represent  numeric values. Let's say your existing color function is C, and your coordinate functions for a parametrized surface are Fx, Fy, Fz.

Gamma:= 1.15:
plot3d(
     [Fx, Fy, Fz],  a..b, c..d, 
     color=  [
          (x,y)-> (1-C(x,y))^Gamma/3, #Hue
          (x,y)-> 1-C(x,y)/4,         #Saturation
          (x,y)-> 1-C(x,y)/7,         #Value
          colortype= HSV
     ],
     lightmodel= NONE,
     style= patchnogrid     
);

There are several parameters that can be adjusted; I've chosen some of them by my personal taste for color .

• Gamma controls the evenness of the distribution between red and green. I gave this one a name because this is a well-known concept (see the Wikipedia article "Gamma correction").
• The 3 in the Hue selects the fraction (1/3 in this case) of the full color spectrum that you want. If you want green to red, it will need to be pretty close to 3.
• The Hue value is subtracted from 1 to make the scale go green to red rather than red to green.
• The 4 in the Saturation controls (to some extent) how "light" the light-green is.
• The 7 in the Value controls (to some extent) how dark the dark-red is (lower values will make it darker).
• lightmodel= NONE is used so that the colors will not change due to shadows when the plot is rotated.

Here is the resulting sub-spectrum:

Brian wrote:

the above solution violates my assumption phi should be between -1.57 and +1.57. why?

Markiyan's Answer addresses well how to get the results. Regarding why you had the problem: solve generally ignores assumptions (see ?solve,details), sometimes without warning. Instead of using assumptions, you can include inequalities with the equations that you pass to solve, as Markiyan did.

You wrote:

so I dont actually need to asign as I did with xs??

Like virtually all commands in Maple, you do need to assign the results of dsolve in order to use them later.

You wrote:

in ordinary circumstances plot(eval(x(t),t=0..3) would be sufficient?

No, the command is plot(eval(x(t), xs), t= 0..3). It requires the xs. Note that you cannot directly plot the results of a dsolve(..., numeric); for the numeric case, use plots:-odeplot instead.

You wrote:

could I similarly use eval(xs,x(t)) ??

No, the command is eval(x(t), xs). The second argument to this type of eval must be an equation, or a list or set of equations--- which is what xs is in the case of a non-numeric dsolve.

You wrote:

ode3 := diff(y(t), t, t) = y(t)*sqrt((1-y(t))^2-1)/(1-y(t))

ics3 := y(0) = -2, (D(y))(0) = 0

ys := dsolve({ics3, ode3}, numeric)

I can ode plot this as the first one.

But my next step was to parametrically plot both solutions like

plots:-odeplot([xs(t), ys(t), t = 0 .. 3])

plots:-odeplot([xs, ys, t = 0 .. 3])

plots:-odeplot([x(t), y(t), t = 0 .. 3])

Unfortunately, each invocation of plots:-odeplot can only be used with the results of a single dsolve(..., numeric). But there are two ways around this problem. The first is to apply dsolve to the ODEs together as a system:

Sol:= dsolve({ode3, ics3, ode, ics}, numeric):
plots:-odeplot(Sol, [[t,x(t)], [t,y(t)]], t= 0..3);

The second is to combine two plots with plots:-display (which can combine any number of plots):

plots:-display([
     plots:-odeplot(xs, t= 0..3),
     plots:-odeplot(ys, t= 0..3)
]);

The second method gives you more flexibility.

You wrote:

how can one plot:

plots:-odeplot(f(xs), t= 0..3);

with some function f

Like this:

plots:-odeplot(xs, [t, f(x(t))], t= 0..3);

The command that you need is

B:= copy(A);

See ?copy for some details why. A few types of "container" structures in Maple are considered "mutable": tables, rtables (which includes Arrays, Matrices, and Vectors), procedures, and modules (which includes Records). Mutable objects need to be fully copied in order for an assignment to one to not affect the other. Note that most structures in Maple are not mutable even if they are containers: lists, sets, sequences.

As you have written it, it's an ODE, not a PDE, because all the derivatives are taken with respect to (wrt) x. Perhaps you meant for the derivative on the right side of the equals sign to be wrt t instead? Since it's a second-order ODE wrt x, two boundary conditions are enough, hence the error message.

The result returned by your dsolve, xs, is of the form x(t) = ....  The plot command is complaining because you passed it an equation. The actual solution is the right-hand side of the equation. There are many ways  to isolate the right-habd side. One way that's easy to remember is rhs(xs). The way I generally prefer is eval(x(t), xs). Under ordinary circumstances you could give the plot command

plot(eval(x(t), xs), t= 0..3);

But these are not ordinary circumstances. The solution, although it appears short, is so complicated that Maple is having a lot of trouble evaluating it numerically. What we need to do is get a numeric solution from dsolve and then use plots:-odeplot to plot it:

xs:= dsolve({ics, ode}, numeric);
plots:-odeplot(xs, t= 0..3);

By the way, the right side of your ode has the useless coefficient 1^2. Did you intend for that to be something else?

How about this for the graphic?

N:= 96:
SP:= [seq](
     plots:-spacecurve(
          [1, k*3*Pi/N+Pi/2/N*floor(k/N), phi], phi= -Pi..Pi, coords= spherical,
          thickness= 9, 
          transparency= 1-(N+k)/3/N
     ), k= 0..N
):
plots:-display(
     [plots:-display(
          [seq](plots:-display(SP[N+2-k..N+1]), k= 1..N+1)
         ,insequence 
     ),
     plots:-tubeplot([0,0,t], t= -1..2, color= blue, radius= 0.05, transparency= .75),
     plots:-tubeplot([t,0,-1], t= -2..2, color= black, radius= 0.05, transparency= .6),
     plots:-tubeplot([0,t,-1], t= -2..2, color= black, radius= 0.05, transparency= .6),
     plot3d(-1, x= -2..2, y= -2..2, color= black, transparency= .5)
     ], axes= none, scaling= constrained, lightmodel= light4, style= patchnogrid,
     glossiness= 1
);

In the future, use square brackets for indexing (such as a[3]), not parentheses (such as a(3)).

Sys:= {seq}(eq(k), k= 1..18):

# Change to square brackets:
Sys:= subs([seq](a(k)= a[k], k= 1..18), Sys):
fsolve(Sys, {seq}(a[k], k= 1..18));

{a[1] = 0.000126988919895945, a[2] = -0.00105660821302691,
 a[3] = 0.00308434899631239, a[4] = -0.00423918949116484,
 a[5] = 0.00282464668788845, a[6] = -0.000740204565874612,
 a[7] = 0.290738925426173, a[8] = -1.17503115010852,
a[9] = 1.66075257813012, a[10] = -0.945945157740450,
 a[11] = 0.174368773622886, a[12] = -0.000158583197276566,
 a[13] = -0.00000150252845025158, a[14] = -0.0968979999980298,
 a[15] = 0.293699162347833, a[16] = -0.332030736202987,
 a[17] = 0.157521639797399, a[18] = -0.0248287447073636}

That took 3.45 seconds on my computer.

The residuals are not very good. You should use a high value of Digits (like 30) for the computation, and you should use more digits in your input constants.  I got good residuals (~ 10^(-19)) at Digits = 30.

This is difficult to diagnose without seeing the output of the previous commands. Nonetheless, I have an idea. Change the line

DE_C_tank_1 := diff(C_1(t), t)+m_out_tank_1*C_1(t)/M_tank_1-C_lost_tank_1 = Charge_rate_V3*C_V3/M_tank_1+Charge_rate_V4*C_V4/M_tank_1;

to

DE_C_tank_1 := diff(C_1(t), t)+~m_out_tank_1*C_1(t)/M_tank_1-C_lost_tank_1 =~ Charge_rate_V3*C_V3/M_tank_1+Charge_rate_V4*C_V4/M_tank_1;

The only difference is that I added two tildes (~), one after the first plus sign and one after the second equals sign. The purpose of the tilde is to "map" the operation over the whole Array. So, the output of the corrected command will be a single Array.

Let me know how that goes.

How about this?

Student:-Calculus1:-VolumeOfRevolution(
     sqrt(1-rho^2), -sqrt(1-rho^2), rho= 0..1,
     axis= vertical,
     output= animation
);
Student:-Calculus1:-VolumeOfRevolution(
     sqrt(a^2-rho^2), -sqrt(a^2-rho^2), rho= 0..a,
     axis= vertical,
     output= integral
);
value(%) assuming a>0;