It's easy to do without calling  LinearAlgebra  package like in the example below:

A:=<1,2; 3,4>; # Matrix of the system 
x:=<x1, x2>; # Vector of unknowns
b:=<5,6>; # Vector of constant terms

System:=convert(A.x=~b, list);  # The system

Presentation in the traditional form:

%piecewise(seq(op([``,System[i]]), i=1..nops(System)));
                    

Edit.

a := x+2*y = 3:
plot(solve(a, y), x = -5 .. 5, -5 .. 5, style = point, symbol=solidcircle, symbolsize=15, numpoints=10, adaptive=false);
                 

You can solve your task by analogy with the following simpler examples:

plots:-implicitplot([seq(y=a*x^2, a=[1,2,3])], x=-1..1, y=0..3, color=[red,blue,green], scaling=constrained);
plot([seq([t, a*t^2, t=-1..1], a=[1,2,3])], color=[red,blue,green], scaling=constrained);

Pay attention to  seq  command. This is one of the most useful and often used commands in Maple.
 

plots:-implicitplot3d(x=y, x= 0..0.01, y= 0..0.01, z=0..100, style=surface, color=grey, transparency= 0.4); 
# Or
plots:-display(plottools:-polygon([[0,0,0],[0,0,100],[0.01,0.01,100],[0.01,0.01,0]], color=grey, transparency= 0.4));
 

Using the symmetry of the ball, the problem reduces to counting the number of points with integer coordinates in the first octant (the result is multiplied by 8).
The procedure named  Onion  returns 2 numbers: the first number is the total number of pieces, the second number is the number of cubes.

The code:

Onion:=proc(R::positive)
local R0, N, n, x, y, z;
R0:=`if`(R::integer,R-1,floor(R));
N:=0; n:=0;
for x from 0 to R0 do
for y from 0 to R0 do
for z from 0 to R0 do
if x^2+y^2+z^2<R^2 then N:=N+1;
if (x+1)^2+(y+1)^2+(z+1)^2<=R^2 then n:=n+1;
fi; fi;
od: od: od:
8*N, 8*n;
end proc:

Examples of use:

Onion(2);
Onion(2.5);
Onion(3);
Onion(10);
Onion(100);
                                 

                                 

1. Your code is correct and answers the point (b). Just add the option  view=[-5..5, -0.5..2.5]  to your code for more visibility.

2. For (a) just replace the option  output=plot  by  output=line  in your code .

3. The code of the animation for  (c) :

restart;
with(Student[Calculus1]):
T:=proc(x0)
local A, B, C, f, k;
uses plots, plottools;
f:=x->2/(1+exp(-x));
k:=D(f)(x0);
A:=line([x0-5/sqrt(k^2+1),f(x0)-5*k/sqrt(k^2+1)], [x0+5/sqrt(k^2+1),f(x0)+5*k/sqrt(k^2+1)], color=red);
B:=plot(f(x), x=-10..10,-2..3.5, color=blue);
C:=disk([x0,f(x0)], 0.08, color=red);
display(A,B,C, scaling=constrained);
end proc:

plots:-animate(T,[x], x=-5..5, frames=90, size=[1000,300], paraminfo=false);

Output:

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);