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

We use the classical method of finding critical points and their subsequent investigation with the help of the second derivative. Maple returns the results in thу range  -Pi .. Pi . If you want the results to be in the range  0..2*Pi , just add  2*Pi  to the corresponding variable:  

restart;
f:=(x,y)->1+8*cos(1/2*x-1/2*y)*cos(1/2*x)*cos(1/2*y);
[solve({D[1](f)(x,y)=0, D[2](f)(x,y)=0},explicit)];
S:=[seq([Student:-MultivariateCalculus:-SecondDerivativeTest(f(x,y),[x,y]=eval([x,y],%[i]))], i=1..nops(%))];
select(s->s[1]<>(LocalMin = []), S);

The final results:

with(plots):
K:=9;
deG:=diff(theta(t),t,t) + mu*diff(theta(t),t)+K*sin(theta(t))= 0;
deL:=diff(theta(t),t,t) + mu*diff(theta(t),t)+K*theta(t)= 0;
Iv:=theta(0)=0.75, D(theta)(0)=2.0;
dom1:=t=0..10;
soln1a:=dsolve({eval(deL,mu=0),Iv});
gr1a:=plot(eval(theta(t), soln1a), dom1, color=blue);
 

Your system, in addition to the unknown, contains several parameters. You can obtain explicit solutions depending on the values of these parameters as follows:

Download System_with_parameters.mw

x^2:
f:=unapply(%, x);
f(t);
f(3);

We can find all 3 solutions as follows:

restart;
z:=x+I*y:
evalc(abs(z)*(z-4-I)+2*I = (5-I)*z);
solve({Re(lhs(%))=Re(rhs(%)), Im(lhs(%))=Im(rhs(%))}) assuming real;
evalf(%);
                     

Edit. I do not know why  solve  (as OP did) fails with this example. Probably this should be seen as a bug.