Yes, we can plot the entire graph without using differential equations according to the following plan:

1. First, we plot the graphs of  f , g , h of the implicit functions, from which we find the maximum range for  t - variable, so that there are real values of  x , y , z.

2. Then, with step = 0.01, we solve the equations for these variables. For every  t  we get  2 values of  x , 2 values of  y  and 1 value of  z .

3. Thus we get 4 branches of the entire graph and plot them.

The final result:

All the codes in the attached file.

curves_4.mw

x , y, z  variables are specified by the equations f , g , h  implicitly. These equations have real solutions not for every  t = 0 .. 2*Pi. This is clearly seen from the graphs of implicit functions. For example, from the graph  f , we see that for  t>3.48...  and  t<2.2...  there are no real solutions:

plots:-implicitplot(f, t=0..2*Pi, x=-10..10, gridrefine=5, view=[0..2*Pi, -1.5..1.5]);
                         

Edit.

You can get an symbolic approximate value of this integral by expanding the integrand in a series in powers of  x :

restart;
S:=12*x^3*c[2]+6*x^2*c[1]+x^2*exp(x^3*c[2])*exp(x^2*c[1]):
P:=convert(series(S, x=0, 10), polynom);
int(P, x = 0..1);
                             

In  s:=LowerCase(s):  line of your procedure there are 2 errors at once:

1. You can not assign anything to formal parameters in the body of a procedure.

2. Strings can not act as names.
 

Here are 2 examples:
It does not work:
P:=proc(x::string)
uses StringTools;
x:=UpperCase(x);
%;
end:
P("y");

But it works:
P1:=proc(x::string)
local s;
uses StringTools;
s:=UpperCase(x);
s;
end:
P1("y");

Edit.
 

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