In addition to the plotting of the ellipse was asked to find a canonical equation. The simplest way to do this - using of the formulas from here

restart;

local D;

f:=(x,y)->3*x^2-3*x*y+6*y^2-6*x+7*y-9;

coeffs(f(x,y));

A,B,C,D,E,F:=%;   # Assigning names to coefficients

theta:=1/2*arctan(B/(A-C));   # Finding the angle of rotation of the canonical system of coordinates relative to the original  coordinate system

solve({-2*A*xc-B*yc=D, -B*xc-2*C*yc=E});   # Finding of the center of the ellipse

assign(%);   # Assigning names to coordinates of the center

# The connection between the canonical  coordinates  and the original  coordinates

x:=xcan*cos(theta)-ycan*sin(theta)+xc;   

y:=xcan*sin(theta)+ycan*cos(theta)+yc;

Eq:=simplify(expand(f(x,y)));   # The equation of the ellipse in canonical coordinates

xcan^2/simplify(sqrt(1/coeff(Eq,xcan^2)*(-tcoeff(Eq))))^`2`+ycan^2/simplify(sqrt(1/coeff(Eq,ycan^2)*(-tcoeff(Eq))))^`2`=1;  # Canonical equation

Visualization:

xcan:=plot(yc+tan(theta)*('x'-xc), 'x'=-2..3.5, color=black):

ycan:=plot(yc-1/tan(theta)*('x'-xc), 'x'=0.1..1.5, color=black):

Ellipse:=plots[implicitplot](f('x','y'),'x'=-2..3.5,'y'=-2..1.5,color=red,thickness=2, gridrefine=5):

labels:=plots[textplot]([[0.4,1.3,"ycan"],[3.2,0.75,"xcan"]], font=[TIMES,ROMAN,14]):

plots[display](xcan,ycan,Ellipse,labels,scaling=constrained);

plot(f(t,a0,b0,c0),  t=0..10);  # or

plot(eval(f(t,a,b,c), {a=a0, b=b0,c =c0), t=0..10);

restart;

Puis:=proc(X, n::integer)

option remember;

if n=0 then 1 elif

n>0 then Puis(X, n):=Puis(X, n-1)*X  elif

n<0 then Puis(X, n):=Puis(X, n+1)/X  fi;

end proc:

Examples of use:

Puis(4,0),  Puis(4,3),  Puis(4,-3),  Puis(4,6),  Puis(4,-6);

                           1, 64, 1/64, 4096, 1/4096

Addition: I understood the meaning of the problem is not to use the sign  ^

rhs  is superfluous. Should be

a := .1994;

modfit3 := a*x^1.5;

f := unapply(modfit3, x);

Compare:

rhs(b);

   Error, invalid input: rhs received b, which is not valid for its 1st argument, expr

A := plot3d([2*cos(phi), 2*sin(phi), z], z = 0 .. 5, phi = 0 .. 2*Pi):

B := plot3d([8*cos(phi), 8*sin(phi), z], z = 10 .. 20, phi = 0 .. 2*Pi):

C := plot3d([(2+6*((z-5)*(1/5)))*cos(phi), (2+6*((z-5)*(1/5)))*sin(phi), z], z = 5 .. 10, phi = 0 .. 2*Pi):

plots[display](A, B, C, axes = normal, scaling = constrained);

Edited. For some reason, some of the code disappears. Now is all right. 

In the line  phiB1 = ...   should be   phiB1 := ...

In fact, your function  omega  depends only on  x  and  the function  y  so

omega:=(x, y)-> (y(x)*(diff(y(x), x))/x);

Example:

omega(x, cos);

                              -cos(x)*sin(x)/x

You can just assign any valid names to specific  styles. Also multiple plots should be given a list and not a set (in order to options conform to plots).

Example:

`1` := dash: `2` := dashdot:

plot([x^2, x+1], x = -2 .. 2, linestyle = [`1`, `2`]);

You must specify the rule that you are going to compare the real and imaginary numbers.

Example:

sort([1+3*I, 2-I, 3], (a,b)->is(abs(a)<=abs(b)));

                        [2-I, 3, 1+3*I]

Read help on  sort  command.

restart;

Lines:=[y, x, x+y-3, x-y+7]:

Points:={[1,1],[4,3],[1,11],[-2,9],[-1,5],[-2,2],[-8,4],[-9,-1],[-3,-3],[3,-3],[6,-1]}:

Points1:=map(t->op([[t[1]-0.5,t[2]],[t[1]+0.5,t[2]]]), Points):

Packs:=combinat[choose](Points1,4):

k:=0:

for p  in Packs do

if convert([seq(add(sign(eval(Lines[i],[x=p[j,1],y=p[j,2]])),j=1..4)=0,i=1..4)],`and`) then

k:=k+1; Teams[k]:=p fi;

od:

Teams:=convert(Teams,list):

L:=ListTools[Categorize]((t1,t2)->convert([seq([seq(sign(eval(Lines[i],[x=t1[j,1],y=t1[j,2]])),j=1..4)]=[seq(sign(eval(Lines[i],[x=t2[j,1],y=t2[j,2]])),j=1..4)],i=1..4)],`and`),Teams):

n:=nops([L]);

Confs:=[seq(L[i,1], i=1..n)]; 

                                                                             n := 13

Confs := [{[-9.5, -1], [-8.5, -1], [3.5, 3], [4.5, 3]}, {[-9.5, -1], [-8.5, 4], [3.5, 3], [5.5, -1]}, {[-9.5, -1], [-3.5, -3], [.5, 11], [3.5, 3]}, {[-9.5, -1], [-2.5, 2], [.5, 11], [5.5, -1]}, {[-9.5, -1], [-2.5, 9], [.5, 1], [5.5, -1]}, {[-9.5, -1], [-2.5, 9], [2.5, -3], [3.5, 3]}, {[-9.5, -1], [-1.5, 5], [.5, 11], [2.5, -3]}, {[-8.5, 4], [-7.5, 4], [5.5, -1], [6.5, -1]}, {[-8.5, 4], [-3.5, -3], [.5, 11], [5.5, -1]}, {[-8.5, 4], [-2.5, 9], [2.5, -3], [5.5, -1]}, {[-3.5, -3], [-2.5, -3], [.5, 11], [1.5, 11]}, {[-3.5, -3], [-2.5, 9], [.5, 11], [2.5, -3]}, {[-2.5, 9], [-1.5, 9], [2.5, -3], [3.5, -3]}]

Visualization:

A:=plots[implicitplot](Lines,x=-10..7,y=-5..12, thickness=3, color=blue):

B:=[seq([seq(plottools[disk](Confs[i,j],0.2,color=red), j=1..4)],i=1..13)]:

plots[display](Matrix(3,4,[seq(plots[display](A,op(B[i]), axes=none), i=1..12)]));

plots[display](A,B[13], axes=none);

Lines:=[y, x, x+y-3, x-y+7]:

Points:=[[1,1],[4,3],[1,11],[-2,9],[-1,5],[-2,2],[-8,4],[-9,-1],[-3,-3],[3,-3],[5,-1]]:

Packs:=combinat[choose](Points,4):

k:=0:

for p  in Packs do

if convert([seq(add(sign(eval(Lines[i],[x=p[j,1],y=p[j,2]])),j=1..4)=0,i=1..4)],`and`) then

k:=k+1; Teams[k]:=p fi;

od:

Teams:=convert(Teams,list);  nops(%);

Teams := [[[1, 11], [4, 3], [-3, -3], [-9, -1]], [[1, 11], [-2, 9], [-3, -3], [3, -3]], [[1, 11], [-8, 4], [-3, -3], [5, -1]], [[1, 11], [3, -3], [-1, 5], [-9, -1]], [[1, 11], [-9, -1], [5, -1], [-2, 2]], [[4, 3], [-2, 9], [3, -3], [-9, -1]], [[4, 3], [-8, 4], [-9, -1], [5, -1]], [[-2, 9], [-8, 4], [3, -3], [5, -1]], [[-2, 9], [1, 1], [-9, -1], [5, -1]]]

                                                                                      9

Visualization of all the 9 solutions:

A:=plots[implicitplot](Lines,x=-10..6,y=-4..12, thickness=3, color=blue):

B:=[seq([seq(plottools[disk](Teams[i,j],0.2,color=red), j=1..4)],i=1..9)]:

plots[display](Matrix(3,[seq(plots[display](A,op(B[i]), axes=none), i=1..9)]));

One equation  a = sin(w*t1 + phi)  is equivalent to an infinite number of linear equations  w*t1 + ph=(-1)^n*arcsin(a)+Pi*n , n is an integer. I take only n=0

The exact solution of the first two equations:

restart;

solve({w*t1 + phi=arcsin(a), w*t2 + phi=arcsin(b)}, {w,phi});

                   {phi = -(arcsin(a)*t2-t1*arcsin(b))/(t1-t2), w = 1/(t1-t2)*(arcsin(a)-arcsin(b))}

If the data obtained as a result of the experiment, the overdetermined system can be solved by the method of least squares.  Maple finds a solution for which the sum of the squares of the differences of the left and right parts is minimal.

LinearAlgebra[LeastSquares]({w*t1 + phi=arcsin(a), w*t2 + phi=arcsin(b), w*t3 + phi=arcsin(c)}, {w,phi}) assuming t1::real,t2::real,t3::real;

{phi = -1/2*(arcsin(a)*t1*t2+arcsin(a)*t1*t3-arcsin(a)*t2^2-arcsin(a)*t3^2-arcsin(b)*t1^2+arcsin(b)*t1*t2+arcsin(b)*t2*t3-arcsin(b)*t3^2-arcsin(c)*t1^2+arcsin(c)*t1*t3-arcsin(c)*t2^2+arcsin(c)*t2*t3)/(t1^2-t1*t2-t1*t3+t2^2-t2*t3+t3^2), w = 1/2*(2*t1*arcsin(a)-arcsin(a)*t2-t3*arcsin(a)-t1*arcsin(b)+2*t2*arcsin(b)-t3*arcsin(b)-t1*arcsin(c)-t2*arcsin(c)+2*t3*arcsin(c))/(t1^2-t1*t2-t1*t3+t2^2-t2*t3+t3^2)}

otat := proc(a::`=`, b::list)

     local q;

     q := isolate(a, y);

seq(eval(q, b[i]), i=1..nops(b));   # returns values of  y  for all values of  x  in the list  b    

end proc:

Example:

otat(x-y = 1, [x = 10, x = 11, x=14]);

                  y = 9, y = 10, y = 13

Addition:  The line  seq(eval(q, b[i]), i=1..nops(b));   can be written shorter  eval~(q, b)[];   with the same output

Boundary problems can rarely be solved symbolically. And your task even has  a parameter K. Such problems are solved numerically for a particular value of the parameter:

K:=1:

Eq := diff(y(x), x, x) = -(x^2+1)*y(x)+K;

sol:=dsolve({Eq, y(-1) = 0, y(1) = 0}, y(x),numeric);

plots[odeplot](sol,[x,y(x)], x=-1..1, color=red, thickness=2);

Here are another 6 ways (from the shortest to the longest) for the same list  L:=[7, 9, 13, evalf(Pi), a, b]:

L^~2;

`^`~(L,2);

map(`^`,L,2);

[seq(t^2, t=L)];

[L[i]^2 $ i=1..nops(L)];

[seq(L[i]^2, i=1..nops(L))];