What 2 lines did you mean? In the plot from your worksheet, I see only one straight line and 2 coordinate axes.
General advice: to find the intersection point of two lines  F(x,y) = 0  and  G(x,y) = 0 , solve the system

solve({ F(x,y) = 0,  G(x,y) = 0}, {x,y});

Edit. Sorry, I didn't notice your second plot right away. So do as Carl wrote or just

restart;
solve({y=5*x-1, y=-x+7});

restart; 
P := 2^n-1: 
select(isprime, [seq(eval(P, n = k), k = 1 .. 100)]);

  [3, 7, 31, 127, 8191, 131071, 524287, 2147483647, 2305843009213693951, 618970019642690137449562111]

Edit. If you do not need prime numbers of the form  P  themselves, but exponents  n , then below is another code:

restart; 
P := 2^n-1: 
select(p->isprime(eval(P, n=p)), [$1..100]);

                 [2, 3, 5, 7, 13, 17, 19, 31, 61, 89]

restart;
a:=3: R:=7: z:=sqrt(R^2-(r-a+R)^2):
plot3d([[r,phi,z],[r,phi,-z]], phi=0..2*Pi, r=2..3, coords=cylindrical, scaling=constrained);

The formula  sqrt(R^2-(r-a+R)^2)  was obtained by elimination  phi  from the system  { z=R*cos(phi), r=a-R*(1-sin(phi)) } 

Edit. I added the option  scaling=constrained so that you can watch the form change depending on the parameters change.

In my opinion the simplest way is to use a triple for-loop:

restart;
C := y^2*z + y*z^2 = x^2:
for x in {0,1} do
for y in {0,1} do
for z in {0,1} do
if C then print(['x','y','z']=[x,y,z]) fi;
od: od: od:

To plot a 3-variable equation, you must use  plots:-implicitplot3d  command. To get the best quality we can first express  z  through  x  and  y and then plot:

restart;
eq:= x^2-y^2*z-y*z^2;
plots:-implicitplot3d(eq,x=0..1,y=0..1,z=-7..7, style=surface, grid=[40,40,40]);
solve(eq, z);
plot3d([%], x=0..1,y=0..1, style=surface);

Edit. A more interesting picture is obtained if we increase the ranges for  x  and  y :  x=-1..1, y=-1..1. It is worth noting one drawback of  plot3d  command in this example: Maple does not draw a surface near z-axis. This is due to the fact that for small values of  y , the value of  z  tends to infinity. The use of  plots:-implicitplot3d  command does not have this drawback. We get good quality when used together  plot3d  and  plots:-implicitplot3d :

restart;
eq:= x^2-y^2*z-y*z^2;
A:=plots:-implicitplot3d(eq,x=-1..1,y=-1..0.01,z=-5..5, style=surface, grid=[50,100,50]):
Sol:=solve(eq, z);
B:=plot3d([Sol], x=-1..1,y=-1..1, view=-5..5, style=surface, axes=normal, grid=[500,500]):
plots:-display(A,B, view=[-1..1,-1..1.2,-5..5]);

Let's help Maple a little. If  t^2 - 4 = s^2 - 4  and  t<>s  then  s = -t . Below we get that there will be infinitely many double points:

restart;
x:=t->sin(t): y:=t->t^2-4:
solve({sin(t)=sin(s),s=-t}, AllSolutions);
eval(%, _Z1=k);
[x,y]=eval([x(t),y(t)], %) assuming k::posint;  # All the double points
# k = 1, 2, 3, ...

                               [x, y] = [0, Pi^2*k^2-4]

The numerical solution:

restart;
x:=t->sin(t): y:=t->t^2-4:
Sys:=[x(t)=x(s),y(t)=y(s)]:
plots:-implicitplot(Sys, t=-10..10,s=-10..10, color=[red,blue], gridrefine=3);
fsolve(Sys, {t=-4..-1,s=1..4});
fnormal(eval([x(t),y(t)], %));

                                    [0., 5.869604404]

Because there are many solutions, then the plot is needed to isolate the root.

With local, you made a new name (it's not clear why?) from the command  pdetest . Just replace the line local pdetest := (ans, de1)  with the line  test := pdetest(ans, de1)

restart;
did:=D^14+23*D^13+144*D^12-30*D^11;
subsindets(did,symbol^integer,t->(D@@degree(t))(y)(x));
convert(%, diff);

You can use a common  if ... fi:  for this:

restart;
str:="A";
x:=10;
if x=10 then str:=cat(str," it was 10"); x:=11 else
str:=cat(str," it was not 10"); x:=9 fi;

Edit.

See help on the  evala  command. Also the  Algebraic  package may be useful for you.

You forgot to call the package  plots . Insert  with(plots):  at the beginning of your worksheet.

Replace in the code the dot  .  with the multiplication sign  * .

The corrected file 

test1_03_new.mw

The error you receive says that Maple does not have the  FourierSeries  package.

You can freely download the  OrthogonalExpansions  package from the Maple Application Center. This package contains a number of commands for working with Fourier series.

https://www.maplesoft.com/applications/view.aspx?SID=33406

(rhs@op)~([rootsq0]);

There will be infinitely many positive solutions too. Below we find the smallest:

restart;
eq := 61*x^2 + 1 = y^2;
[isolve(eq)];
map(t->eval([x,y],t),%);
map(t->expand(eval(t,_Z1=1)), %);
[x,y]=~select(t->t[1]>0 and t[2]>0, %)[];

                     [x = 226153980, y = 1766319049]

Edit. It is convenient to find the desired solutions using the procedure:

restart;
Eq := 61*x^2 + 1 = y^2:
Sol:=n->subsindets(eval(isolve(Eq)[4],_Z1=n),realcons,expand):

# Examples of use:
Sol(1);
Sol(2);