f:=x->x^2*sqrt(25-x^2):
a:=0.: b:=5.: n:=100: h:=(b-a)/n: S:=0:
for j from 0 to n-1 do
S:=S+(f(a+j*h)+f(a+(j+1)*h))/2*h;
end do:
S;

restart;
with(plots):
local gamma:
g:=2:
rr:=3:
V1:=arrow([0.5,0], [-4.5,0], width=[0.006, relative=false], head_width=[0.07, relative=false], head_length=[0.2, relative=false], view=[-4..4,-1..1]):
V2:=arrow([1.5,0],[2.5,0], width=[0.006, relative=false], head_width=[0.07, relative=false], head_length=[0.2, relative=false], view=[-4..4,-1..1]):
Seg:=plot(0, x=1/g+0.1..rr/g-0.1, color=black, thickness=7):
C:=plot([[0.55,0],[1.45,0]], style=point, color=black, symbol=circle, symbolsize=14):
L:=plot([[0,t,t=-0.04..0.04],[1,t,t=-0.04..0.04],[2,t,t=-0.04..0.04]], color=black, thickness=5):
T:=textplot([[0,-0.15,0],[1,-0.15,1],[0.55,-0.15,1/gamma],[1.45,-0.15,r/gamma],[2,-0.15,r]], font=[times,bold,16]):
display(V1,V2,Seg, C, L, T, size=[800,400], axes=none);

   Output:

restart;
y1:=-1;
y2:=eval(y, solve({y=x, x=2-y^2, x>0}));
int(1, [x=y..2-y^2, y=y1..y2]);
                                                            

x,y := 1,2;
x%+y=x+y;
                                               

Replace  sum  by  add . Now it works:

restart;                                                                  
with(plots):                                                               
g:=(x,y)->sin(Pi*(x+y))*(x^2+y^2);                                        
q:=(n,k,l,x,y)->binomial(n,k)*binomial(n-k,l)*(((1+x)/2)^k)*(((1+y)/2)^l)*(1-((1+x)/2)-((1+y)/2))^(n-k-l) ;

 K:=(n,x,y)->(((1+n)*(2+n))/4)*add(add(q(n,k,l,x,y)*int(int(q(n,k,l,t,u)*g(t,u),u=-1..-t),t=-1..1),l=0..n-k),k=0..n);  

p1:=plot3d(g(x,y),x=-1..1,y=-1..-x,color=blue):                           
p2:=plot3d(K(1,x,y),x=-1..1,y=-1..-x,color=brown):                        
p3:=plot3d(K(2,x,y),x=-1..1,y=-1..-x,color=yellow):                       
p4:=plot3d(K(5,x,y),x=-1..1,y=-1..-x,color=green):                        
p5:=plot3d(K(10,x,y),x=-1..1,y=-1..-x,color=red):                                                              
display([p1,p2,p3,p4,p5]);
 

restart;
int(eval(2*x^i*(x+n)^m*sqrt(-x^2+x), n=0), x = 0 .. 1);
int(2*x^i*(x+n)^m*sqrt(-x^2+x), x = 0 .. 1) assuming n>=1;
                                      

f := x->exp(x^5):  n := 3: 
f(0)+convert(f(x)-f(0), FormalPowerSeries):
evalindets(%, 'factorial(anything)', t->``(op(1,t))!):
subs(infinity = n-1, %): 
value(%);

                               

Edit.

restart; with(plots):
H := 11;
B := x; W := x;
t__E := (13100*(B*B)*H*(.23-.10)*1.2)/(1500*alpha);
A__t := 2*B*W+2*H*(B+W);
OF := alpha/A__t;
for alpha from 0.1 by 0.1 to 1.8 do XX[alpha] := solve({0 < x, OF < 0.6951871658e-2, t__E < 200}, x); XY1[alpha]:=[alpha,op([1,1],XX[alpha])];
XY2[alpha]:=[alpha,op([2,2],XX[alpha])];
 end do;
XY1:=sort(convert(XY1,list), (a,b)->a[1]<=b[1]); XY2:=sort(convert(XY2,list), (a,b)->a[1]<=b[1]);
plot([XY1, XY2], color=[blue,red]);

Because you did not submit the text of the code (only pictures), I show how this can be done with the example of 3 spheres:

g1:=x^2+y^2+z^2-16:  g2:=(x+1)^2+(y+3)^2+z^2-16:  g3:=x^2+y^2+(z+2)^2-16:
plots:-implicitplot3d(max(g1,g2,g3)=0,  x=-4.5..3.5, y=-4.5..2.5,z=-4.5..4.5, style=surface, axes=normal, scaling=constrained, lightmodel=light4, numpoints=1000000); 
 

Edit.

Did you forget to call  NumberTheory  package? See help on  NumberTheory:-ContinuedFraction  command. The parameters  q  and  t  should be specified.

restart;
y:=x->5*x^2:  x0:=-1:  # Specify the function and the point
Dy:=limit((y(x0+t)-y(x0))/t, t=0);  # Find the derivative at the point x0
                                         

is(y(x0)=5);  # Check that the point (-1,5) is on the graph of y=5*x^2 
Tang:=y(x0)+Dy*(x-x0);  # Find the tangent line at the point (-1,5)
                                    

plot([y(x), Tang, [[x0,y(x0)]]], x=-2..2, -6..12, color=[red,blue,red], style=[line,line,point], symbol=solidcircle, symbolsize=12, tickmarks=[default,spacing(1)], size=[200,500]);  # The plotting
                     

PS. The derivative was found using the limit (as OP wrote 'Use the formal definition of the derivative..."). Of course, it is shorter to look for the derivative at a point using the differential operator  D : 

D(y)(-1);

Download _The_plotting.mw

In your method, you will always get a graph of the zero function, because your function is nonzero on a discrete set of points. Try another method:

restart;
plot([seq([i,t,t=0..1], i=1..5)], color=blue, view=[0..5.5, 0..1], size=[800,500]);

Here are 3 ways:

a:=2: b:=5:
plot([[0,-a],[0,a],[b,a],[b,-a],[0,-a]], linestyle=2, color=red, thickness=3, scaling=constrained);  # The first
plots:-display(plottools:-curve([[0,-a],[0,a],[b,a],[b,-a],[0,-a]], linestyle=2, color=red, thickness=3), scaling=constrained);  # The second
plots:-display(plottools:-rectangle([0,a],[b,-a], style=line, linestyle=2, color=red, thickness=3), scaling=constrained);  # The third

Edit.

Shorter:

subs([a=2*x,b=3*y], expand(sin(a+b)));

Here are the plots in 2 styles: as points and as a line:

X:=[seq(i, i=0..0.5, 0.05)]; 
Y:=[0.486935382154125, 0.485087274176440, 0.483255914856304, 0.481441076124814, 0.479642533998987, 0.477860068520125, 0.476093463645627, 0.474342507167362, 0.472606990609939, 0.470886709216958, 0.469181461771770];
plot(X,Y, view=0..1, style=point, color=red, symbol=solidcircle, symbolsize=12);
plot(X,Y, view=0..1, color=red, thickness=2);