See help  ?dsolve, numeric, DAE

Addition. Here is an example from the help, in which I added only one line of code for plotting the function  x(t) : 

Download DAE-system.mw

Download Pendulum2time_new.mw

>  restart; r:=1.1241: x0:=-0.1241: y0:=0.0: >  z:=x+I*y: f:=unapply(evalc~([(Re,Im)(z+1/z)]), x,y); F:=plottools:-transform(f): C:=plots:-implicitplot((x-x0)^2+(y-y0)^2=r^2, x=-2..2, y=-2..2, color=red, thickness=2) : plots:-display(<C, F(C)>, scaling=constrained, size=[600,250]);

Addition. The airfoil bend is affected by changing  y0-parameter. The animation below clearly shows this (y0 changes from 0 to 0.1 with the step 0.001):

restart;
r:=1.1241: x0:=-0.1241:
z:=x+I*y:
f:=unapply(evalc~([(Re,Im)(z+1/z)]), x,y);
F:=plottools:-transform(f):
C:=y0->plots:-implicitplot((x-x0)^2+(y-y0)^2=r^2, x=-2..2, y=-2..2, color=red, thickness=2) :
plots:-animate(plots:-display,[ 'F(C(y0))', scaling=constrained], y0=0..0.1, frames=101, size=[800,350]);
               
 

Download Joukowski.mw

restart;
u1,u2,u3,u4 := <1,1,-2,-3>, <-1,-9,b,11>, <-1,3,2,-1>, <b,-10,-4,6>;
A:=<u1 | u2 | u3 | u4>;
solve(LinearAlgebra:-Determinant(A)<>0) ;

You might be interested of the procedure named  IntOverDomain, which  finds a double integral over an arbitrary domain bounded by a non-selfintersecting piecewise smooth curve.  See here

This procedure has a wider scope (applicable to a wider class of regions) than that mentioned in the previous answer.

Your function  Expression  has no critical points, since we have
diff(Expression, a)=80*cos(theta), diff(Expression, b)=30*sin(theta) ,
and the  sin(theta)  and  cos(theta)  are not equal to 0 for the same angle  theta  (because  (sin(theta))^2+(cos(theta))^2=1) .

Therefore, it makes sense only to search for the largest or smallest value of  Expression  on some set. In the example below, we find the largest value of the  Expression  on the set  {a>=0,b>=0,theta>=0,a<=1,b<=1,theta<=1} :

Expression :=30*a*sin(theta)+80*cos(theta)*b;
Optimization:-Maximize(Expression, {a>=0,b>=0,theta>=0,a<=1,b<=1,theta<=1});
      [85.4400374531753073, [a = 1., b = 1., theta = 0.358770670270572]]

Addition. If you look for the maximum only by the changing the variable  theta  (with constants  a  and  b), here Maple makes a serious bug (see this post)

I think that your system with 3 parameters (c, d, k) is advisable to solve numerically for specific parameter values. The procedure  Sol  does this. You can play with it by selecting the desired values of the parameters and getting the corresponding values of  N  and  P .

restart;
Sys:=[c*(N-k)*(1-exp(-P))=P, -d*((k-N)*exp(-P)-k)=N];
Sol:=proc(L::list)
fsolve(eval(Sys,[c,d,k]=~L), {N=0..infinity,P=0..infinity}, avoid={N=0,P=0});
end proc:

#  Example of use:

Params:=[2,3,4]:  # The list of parameters (for example)
Sol(Params); # The solution of the system Sys
plots:-implicitplot(eval(Sys,[c,d,k]=~Params), N=0..2*eval(N,%), P=0..2*eval(P,%), color=[red,blue], thickness=2, gridrefine=3);  # Visualization
               

Download System.mw

For example (for  F  from your worksheet):

F := plot([sigma1, sigma2], sigma = -15 .. 15, color = [RED], thickness = 2, legend = [" &xi;,=2, MCST(1) ", " &xi;=2, MCST(1)"]):
DataF:=plottools:-getdata(F);

# Then

op([1,3], [DataF]);
# Or
op([2,3], [DataF]);
 

See help on  ?background

If now you just want to get the value of the polynomial  v(x)  for the found parameter values  a  and  b, then execute the following in the end of the code:

v(x);
                             (1/2)*x^2-(1/6)*x^3+2/3*x

Your expression  yy  is identically 0 for any  z :

Maple will understand what you want if you do so:

restart;
with(LinearAlgebra):
u := <u1, u2>;
v := <v1, v2>;
A:=<u | v>;
B := Transpose(A).A;

restart;
r := 1.9;
i := floor(r);

You can use  Student:-Calculus1:-Roots  command to find all the real roots in the specified interval. It does not need  maxsols  option and can work symbolically.

Examples:

Student:-Calculus1:-Roots(eq, beta=-10..10);
Student:-Calculus1:-Roots(sin(x), x=-10..10);
 

Of course in the real domain it is better to write  int(tan(x), x) = -ln(abs(cos(x)) . The same for the integrals  int(1/x, x) , int(1/(x^2-1), x)  and so on. But the result of Maple is also correct, since it works by default in the complex domain and the results may differ by a complex constant. See

simplify([ln(-2), ln(-5)]);
                                           [ln(2)+I*Pi, ln(5)+I*Pi]

Since, as you can see, the results have the same complex constant, the final result is correct when calculating certain integrals using the Newton – Leibniz formula:

A:=int(tan(x), x);
B:=int(tan(x), x=2*Pi/3..3*Pi/4);
C:=``(expand(eval(A, x=3*Pi/4))) - ``(eval(A, x=2*Pi/3));
expand(C);
              
                                               

The first result  B  ( cos(x)<0  in this range  2*Pi/3..3*Pi/4 ) is calculated by Maple, and the second one  C  is calculated directly by the Newton-Leibniz formula. The results are the same and correct.

Probably a more complete answer why Maple is designed in this way, the developers of  int  procedure can give you.