One thing to watch: don't use with in a procedure.  If DEtools was implemented as a module, I'd suggest using use instead, but DEtools is still a table-based package. So:

myproc:= proc(v,phi,sirka,k)
  local eq1, eq2, sys; 
  eq1 := diff(x(t), t) = v*cos(phi):
  eq2 := diff(y(t), t) = -v*sin(phi)+k*x(t)*(sirka-x(t)):
  sys := [eq1, eq2]:
  DEtools[DEplot](sys, [x(t), y(t)], t = 0 .. 15, [[x(0) = 0, y(0) = 0]], x = 0 .. sirka, 
      y = -10 .. 10, stepsize = .1, thickness = 1, linecolor = red, arrows = none);
end proc;

For the osculating circle, you might look at RadiusOfCurvature in the Student[VectorCalculus] package.

The osculating plane contains the tangent and normal vectors, and is normal to the binormal vector: see TNBFrame, or TangentVector and PrincipalNormal, in that package.

For the definition of osculating sphere, see mathworld.wolfram.com/OsculatingSphere.html

For example, something like this?

>  with(plots):
   L := [seq([c,rand(10)()],c="a".."m")]; # just some random data
   display([pointplot([seq([L[i,2], nops(L)-i], i=1..nops(L))], symbol=solidcircle, 
              symbolsize=20),
      seq(plot([[0,i],[10,i]],colour=cyan),i=0..nops(L)-1)],
      axes=box, tickmarks=[10, [seq(nops(L)-i=L[i,1],i=1..nops(L))]]);

You might start by looking at polar coordinates.

Unfortunately, Maple does not have good facilities for low-level manipulation of sums.  The sumtools or SumTools package would seem to be good places for these (analogous to IntegrationTools), but these are dedicated to sophisticated methods for finding closed-form summations.  Well, here's a first attempt at writing a procedure for this.  By the way, your rewriting is not quite right: it would require S(0) = 0 (try it for N=1).  The procedure below corrects that.

SumByParts:= proc(E::Sum(algebraic,name=range), S::algebraic)
   local M,q,v,x,a,b,k;
   if op(0,E) = 'sum' then return Parts(Sum(op(E)),S) end if; 
   if not typematch(E, 
       Sum(q::algebraic,x::name=a::algebraic..b::algebraic),'M')
     then error "Wrong type" end if;
   assign(M);
   if type(q,`*`) and member(S,[op(q)]) then v:= q/S
   else v:= normal(q/S) 
   end if;
   Sum((S - eval(S,x=x-1))*sum(eval(v,x=k),k=x..b) 
        + eval(S,x=a-1)*v,x=a..b)
   end proc;
 

For example:

> Q:= Sum((a(j+1)-a(j))*j,j=1..N);
  SumByParts(Q, j);
 

 Sum(a(N+1)-a(j),j = 1 .. N)  

The straightforward way (without using lambda), given n, is to start with the
set of integers {$1..n-1}, select those coprime to n (using select and igcd),
and compute a^(n-1) mod n for all those a.

> Q:=p(x,a) *diff(w(x),x$2) + q(x,a) *diff(w(x),x) + r(x,a)* w(x);
   coeff(Q,diff(w(x),x$2));

> (e,v) := Eigenvectors(subs(cos(x)=c, sin(x)=s, T));
   e:= subs(c=cos(x),s=sin(x),e);

Yes, there should be two solutions, because lhs(f) goes to -infinity as x -> 3/4 but is positive when x is not extremely close to 3/4.  For the equation
1 + x^2 + epsilon*ln(4*x-3) = 0, I think as epsilon -> 0+ the solution is approximately

x = 3/4 + 1/4*exp(-25/16/epsilon)-3/32/epsilon*exp(-25/8/epsilon)

which for epsilon = 1/100 agrees (at Digits = 100) with Alec's second numerical solution.

Speaking of possible changes to if:

The failure of the construction

if a < b then ... end if

when a and b evaluate to real constants that are not of type numeric is a very common cause of grief, and not just to Maple novices.  Always checking conditionals using is instead of evalb might, I suppose, cause performance to suffer.  But using is when Maple would otherwise return the dreaded "cannot determine if this expression is true or false" error shouldn't affect performance when the is is not necessary, and I can't see it breaking much existing code.   So this is one case that I think it would be useful to change.

> with(ExcelTools):
   A:= Import("6629_Density.xls");
  Temperature:= convert(A[4..11,1],list);
  Pressure:= convert(A[3,2..6],list);
  Density:= Matrix(A[4..11,2..6]);
  plots[surfdata]([seq([seq([Temperature[i],Pressure[j],Density[i,j]],j=1..5)],i=1..8)],
      axes=box,labels=['T','P','D']);

> plot(x,x=0..1,legend="foo",legendstyle=[font=[TIMES,ITALIC,20]]);

(works in Standard, not in Classic).

> with(plots):
   eq:=1/2*L*(14.4/(1.5+R))^2<=423e-6;
   implicitplot(eq,L=757.4e-6..1135.1e-6, R = 13.85..15.95, filledregions);

Is that what you wanted?

I think the problem is that, unlike most mathematical objects, Maple can have two identical Vectors in memory, and "=" will not consider them to be the same.
The recommended way to compare them would be with Equal in the LinearAlgebra package.

Why do you need to get the list of points from the plot command?  It would be much simpler to just construct your A directly something like

> A := [seq( [X[i], evalf(f(X[i]))], i = 1 .. nops(X))];

where X is the list of x values you want to use.