plots:-display(plottools:-circle());
plot([cos(t),sin(t), t=0..2*Pi]);
plots:-implicitplot(x^2+y^2=1, x=-1..1,y=-1..1);

More 3 ways:

plot(1, phi=0..2*Pi, coords=polar);
plot([sqrt(1-x^2),-sqrt(1-x^2)], x=-1..1, color=black);
geometry:-draw(geometry:-circle(c,[geometry:-point(A,0,0),1]));

Edit.
 

f:=sqrt(x)*ln(y);
solve(map(Im, evalc(f)), {x,y});
                         

because in accordance with your boundary conditions:

eval(dsys3[3], [w(x) = 0, (D(w))(x) = 0, ((D@@2)(w))(x) = 0]);
                                         0.3693149535

Do not use  and  for this. Use curly braces instead in  solve  command:

solve({x>-infinity, x<infinity, x<>1/2}, x);
                             

You can not plot contour graphs for the function  lambda , because for fixed variables  Nb  and  Gamma2 , you get a function from only one variable  delta2 . Below are contour plots for  lambda  when  delta2=0.02..0.1  and   Nb=0.1...0.3 :

restart;
  h:=z->1-(delta2/2)*(1 + cos(2*(Pi/L1)*(z - d1 - L1))):
  K1:=((4/h(z)^4)-(sin(alpha)/Gamma2)-h(z)^2+Nb*h(z)^4):
  lambda:=unapply(Int(K1,z=0..1), Gamma2):
  L1:=0.2:
  d1:=0.2:
  alpha:=Pi/6:
 with(plots):
  display
  ( Matrix(3,2, [seq(contourplot(lambda(Gamma2), delta2=0.02..0.1,       Nb=0.1...0.3, labels=[typeset(`&delta;1`), typeset(conjugate(`&Delta;p`))], title=typeset("Effect of ", ''alpha'', " when ", Gamma,"2=", Gamma2)), Gamma2 in [10,20,30,40,50,-10])]));  

Purely visually, these graphics are almost indistinguishable. Compare:

contourplot(lambda(-10), delta2=0.02..0.1,  Nb=0.1...0.3);
contourplot(lambda(50), delta2=0.02..0.1,  Nb=0.1...0.3);
 

Since your parameters  N  and  m  each have 3 values, then there will be 3 * 3 = 9 all possible combinations. Therefore, you need nested loops to solve the system. But there is a new problem: Maple writes that there are too few boundary conditions:

Download MHD_cchf_2_new.mw

You need to fill the space between the cube and the sphere with something, but so that it does not block the sphere. I suggest that you do this with red points using  plots:-pointplot3d command. Here are 2 options: in the first variant, the points fill the entire space between the cube and the sphere, in the second one - half of this space:

with(plottools): with(plots):
A:=sphere(color=pink, style=surface):
B1:=pointplot3d([seq(seq(seq(`if`(i^2+j^2+k^2>1,[i,j,k],NULL),i=-1..1,0.1),j=-1..1,0.1),k=-1..1,0.1)], color=red,symbolsize=12):
B2:=pointplot3d([seq(seq(seq(`if`(i^2+j^2+k^2>1  and i<=0,[i,j,k],NULL),i=-1..1,0.1),j=-1..1,0.1),k=-1..1,0.1)], color=red,symbolsize=12):
display(<display(A,B1) | display(A,B2)>);
       

Here are 2 options for plotting: the first 10 and the first 20 Fibonacci numbers. As the numbers increase rather quickly, on the second plot, the difference between the first 9 numbers is almost invisible:

seq(combinat:-fibonacci(n), n=1..20);
plot([seq([n,combinat:-fibonacci(n)], n=1..10)], style=point);
plot([seq([n,combinat:-fibonacci(n)], n=1..20)], style=point);
 

I made a few fixes in your code:

restart;
N:=4; alpha:=5*3.14/180; r:=10; Ha:=5; H:=1;
dsolve(diff(f(x),x,x,x));
Rf:=diff(f[m-1](x),x,x,x)+2*alpha*r*sum*(f[m-1-n](x)*diff(f[n](x),x),n=0..m-1)
+(4-Ha)*(alpha)^2*diff(f[m-1](x),x);
dsolve(diff(f[m](x),x,x,x)-CHI[m]*(diff(f[m-1](x),x,x,x))=h*H*Rf,f[m](x));
f[0](x):=1-x^2;
for m from 1 by 1 to N do
CHI[m]:=`if`(m>1,1,0);
f[m](x):=int(int(int(CHI[m]*(diff(f[m-1](x),x,x,x))+h*H(diff(f[m-1](x),x,x,x))
+2*h*H*alpha*r*(sum(f[m-1-n](x)*(diff(f[n](x),x)),n=0..m-1))+4*h*H*alpha^2*
(diff(f[m-1](x),x))-h*H*alpha^2*(diff(f[m-1](x),x))*Ha,x),x)+_C1*x,x)+_C2*x+_C3;
s1:=evalf(subs(x=0,f[m](x)))=0;
s2:=evalf(subs(x=0,diff(f[m](x),x)))=0;
s3:=evalf(subs(x=1,f[m](x)))=0;
s:={s1,s2,s3}:
f[m](x):=simplify(subs(solve(s,{_C1,_C2,_C3}),f[m](x)));
end do;
f(x):=sum(f[i](x),i=0..N);
hh:=evalf(subs(x=1,diff(f(x),x))):
plot(hh,h=-1.5..-0.2);
A(x):=subs(h=-0.9,f(x));
plot(A(x),x=0..1);

The  coeff  command does not work for extracting the coefficients of polynomials from several variables. A special procedure is required for this. The  coefff  procedure extracts the coefficient in front of a monomial  t of the polynomial  P  from the variables  T .
 

Download coefff.mw

Edit.

In the first equation there should be the multiplication sign after  x , otherwise everything that stands in parentheses (after  x)  disappears.
See:

x:=2:
x(h(t));
                                 2

TM:=proc(N)
if N<2 then error "Should be N>=2" fi;
Matrix(N+1, [[seq(p[i],i=1..N-1),o[3]], [o[1],seq(q[i],i=1..N-1),o[4]], [o[2],seq(r[i],i=1..N-1)]], shape = band[1,1], scan=band[1,1]);
end proc:

Example of use:

TM(5);
                                

To multiply the matrices, I replaced  &*  with the dot  . . Now everything works.

Help_(3)_new.mw

Edit.

for this. 

Example:

plot([cos(t), sin(t), t=-1..1], x=0..2, y=-2..2, color=red, scaling=constrained);   # Or

plot([cos(t), sin(t), t=-1..1], color=red, scaling=constrained, view=[0..2,-2..2]); 
 

plot( x[3]^5, caption = typeset("\n A plot of %1.", x[3]^5), captionfont=[times, 20] );