You can extract data from a plot using the  plottools:-getdata  command. The data is retrieved as a two-column matrix. See a simple example below:

restart;
plot([y^2,y, y=-1..1]);
plottools:-getdata(%);
%[3][1..10];

The integral must be introduced in the inert form:

restart;
A:=Int(exp(-4*tau)*sin(2*tau)/2, tau=0..t);
Student:-Calculus1:-ShowSolution(A);

Download Int.mw

It can be written much shorter using arrow-notation. We see that Maple finds this sum symbolically in a closed form. For large numbers N , this is much more efficient:

restart;
sum2N:=N->sum((N+k)^2, k=0..N);

# Examples
sum2N(N);
sum2N(10);

Alternative :

[seq(rhs~(t)[],t=[a])];

                  [-23, -12, -34, 87, 18, 98, 27, 93, 45, 68]

restart;
Sys := Y(t)=diff(y(t),t,t), diff(y(t), t, t)+y(t)*abs(y(t)) = 0; 
ic1 := y(0) = 1, D(y)(0) = 0; 
dsol1 := dsolve({Sys, ic1}, numeric, range = 0 .. 10);
plots:-odeplot(dsol1,[t,y(t)]); 
plots:-odeplot(dsol1, [t,diff(y(t),t)]);
plots:-odeplot(dsol1, [t,Y(t)]);

The  Picture  procedure is not suitable for your task. Use the built-in maple-procedure  plots:-inequal  for this:

restart;
plots:-inequal({x^2+y^2<=9,x^2+y^2>=4}, x = -4 .. 4, y = -4 .. 4, optionsfeasible = [color = yellow], optionsclosed = [color = brown, thickness = 3]);

In fact, you are trying to find all 3 roots of the third degree of 1 in the complex domain. The shortest way to do this is to simply solve the equation  z^3=1 :

restart;
solve(z^3=1);

That is the  LinearAlgebra:-MatrixExponential  command. See help for details on this command.

If you just type  a^3*sqrt(a);  Maple calculates it as  a^(7/2) . The desired display requires a special procedure. Procedure  P  displays an expression of the type  symbol^fraction  in the desired form.

P:=proc(p::symbol^fraction)
local e, a, n, k, m;
uses InertForm;
e:=op(2,p); a:=op(1,p);
n:=floor(e); k:=e-n; m:=surd(a^numer(k),denom(k));
`if`(n=0,m,`%*`(a^n,m));
Display(%,inert=false);
end proc:

Examples of use:

P(a^(7/2));
P(b^(8/3));
P(a^(2/5));

To avoid this bug, after differentiating with respect to  t  in the resulting series, you can simply write separately the term for k = 1 :

restart;
S:=Sum(epsilon^k*f[k](t)/k!, k=1..infinity);
S1:=diff(S, t);
S1:=eval(op(1,S1),k=1)+Sum(op(1,S1),k=2..infinity);
S2:=simplify(diff(S1,epsilon$2));

The first 5 terms of the series  S2 :

seq(eval(op(1,S2), k=n), n=2..6);

Figure shows a part of the surface inside a sphere of radius 4.5 :

restart;
with(plots):
x0 := x - Pi/4;
y0 := y - Pi/4;
z0 := z - Pi/4;
f := (x,y,z)->10*sin(x0)*sin(y0)*sin(z0) + 10*sin(x0)*cos(y0)*cos(z0) + 10*cos(x0)*sin(y0)*cos(z0) + 10*cos(x0)*cos(y0)*sin(z0) - 7/10*cos(4*x) - 7/10*cos(4*y) - 7/10*cos(4*z);
implicitplot3d(piecewise(sqrt(x^2+y^2+z^2)>4.5, undefined, f(x,y,z)-9),  x=-4.5 .. 4.5, y= -4.5 .. 4.5,  z=-4.5 .. 4.5, grid = [100, 100, 100], style = patchnogrid);

You write your system, first in vector form (as it is defined), and then convert to normal form and solve using  dsolve  command:

restart;
Sys:=diff(<x1(t),x2(t)>,t)=~<-4,-2;6,3>.<x1(t),x2(t)>+<2/(exp(t)+1),-3/(exp(t)+1)>;
dsolve(convert(Sys,set) union {x1(0)=0,x2(0)=ln(2)});

You missed the multiplication sign:

evalf(((12*sqrt(13))/13 - 3/4)/(1 + ((12*sqrt(13))/13)*(3/4)));

                                                 0 .7374399056

Use  plots:-textplot  command for this. For example, add the new line  Anim4  to your code (for the example1):

Anim4 := animate(textplot, [[op(convert(Position(t)+(1.3/4)*Vitesse(t), list)), "V"], align = [right, below], font = [times, bold, 18]], t = t0 .. t1, frames = 100):

When creating complex animations, it is convenient to first write a procedure that creates one frame and then use  plots:-animate  command. Here's the original example:

restart;

P:=proc(t)
local f, g, alpha, beta, Position, Vitesse, Accélération, P1, P2, V1, V2 ;
uses plots, plottools;
f := t->sin(3*t)*cos(t);
g := t->cos(3*t)*sin(t);
alpha := (1/4)*Pi;
beta := 3*Pi*(1/4);
Position := <f(t), g(t)>;
Vitesse := <(D(f))(t), (D(g))(t)>;
Accélération := <(D@@2)(f)(t), (D@@2)(g)(t)>;
P1:=plot([f(s), g(s), s=alpha .. t], color=green, thickness=2);
P2:=plot([[f(t), g(t)]], style=point, symbol=solidcircle, color=blue, symbolsize=20);
V1:=arrow(Position,0.1*Vitesse, width=0.015, color=blue); # Vector of vitesse
V2:=arrow(Position,0.1*Accélération, width=0.015, color=red); # Vector of accélération
display(P1,P2,V1,V2);
end proc:

plots:-animate(P,[t], t=Pi/4..3*Pi/4, frames=90, size=[500,500]);