restart;
with(plottools): with(plots):
s := display(sphere([0,0,0],1), style=surface, transparency=0.8);
c := display(cuboid([-1,-1,-1], [1,1,1]), transparency=0.3, style=surface);
display([c,s], axes=none, orientation=[-65,72,0]);

 

Your calculations assume that the functions { C(n,x), S(m,x) } form an orthogonal family with respect to the weighted integral (weight = 1/x) on the interval (1,2), but that's not true.  Let's check:

We see that the result is not zero in general.

You may already know the difference between the colon (":") and a semicolon (";")  command terminators.  In both cases the command is executed, but in the case of the colon the result is not shown.

Important:  Always use the semicolon terminator when developing new code, especially when you are having trouble with it.  Suppressing the result through a colon deprives you of the assurance that your command works as intended.

With this in mind, go back to your code and replace the colon terminators with semicolons, and then come back here if you still have problems.

Just by looking at your equations it is clear that y(x)≡0 is a solution.  I suppose you are looking for nonzero solutions.

Your equations arise in the analysis of buckling of an Euler beam.  A combination of the parameters P, K(x), and L that yield nonzero solutions correspond to the beam's buckling modes.

You must be familiar with the analysis corresponding to a uniform beam, that is, when K(x) is a constant.  There one shows that buckling modes exist only when L*sqrt(P/K) = nπ, where n is a positive integer.  When K is not a constant, there is no general formula for buckling modes.  The best you can hope for is some numerical experiments for specific choices of K(x).

restart;
m := f__1(x,y)*A + f__2(x,y)*B + f__3(x,y)*C;
f := <coeffs(m, [A,B,C])>;

I don't know how to do this in Maple but that sort of thing is easily done using specialized image manipulation tools.  I would do it with Imagemagick's command-line tools but I am sure there are many other alternatives, perhaps gimp or photoshop.

It takes only two imagemagick commands to produce the desired composite.  Let's call your two images fig.png and fig-zoomed.png.

• Step 1:  Scale fig-zoomed.png  to half its original size and call the result fig-zoomed-scaled.png:
      convert -scale 50% fig-zoom.png fig-zoom-scaled.png
• Step 2: Produce the composite which we call out.png:
      composite -geometry +150+50 fig-zoom-scaled.png fig.png out.png

The "-geometry +150+50" option says that the inset image will be offset relative to the upper-left corner by 150 pixels (horizontally) and 50 pixels (vertically).  Here is the result:

Name the plots, as in
p1 := plot(<a,c,e,g,i,a>, <b,d,f,h,j,b>);
...
p2 := plot(<a,c,e,g,i,a>, <b,d,f,h,j,b>);
and then do:
plots:-display([p1,p2], scaling=constrained);

Download mw.mw

That looks like a bug to me.  Maple 2018 does not crash, but fails with the error message "Error, (in Explore) while processing result".

Here is a workaround that works on Maple 2018.  See if it works for you.  Change
caption = typeset(x^2 + y^2)
to
caption = typeset("%1+%2", x^2, y^2)
 

Download mw.mw

I can't tell by looking in your worksheet what you want your animation to show.

But it appears that you know how to produce the individual frames/graphics of the animation. Let's say these are called p1, p2, ..., pn.  To animate, do:
plots:-display([p1, p2, ..., pn], insequence);

I made some effort to understand the purpose of your calculations and untangle their logic. The result is what I have shown below.  Be sure to check every detail because it is possible that I have misinterpreted your worksheet or introduced typographical errors.

Download mw.mw

Download newton-interpolation.mw

Tomleslie said there are other ways.  Here is one to produce the same object:

restart;
with(plots): with(plottools):
display(cuboid([0,0,0], [1,2,5]), scaling=constrained);

You don't need Maple for that—just looking at the equations it is evident theta(r)≡0 and sigma(r)≡0 is a solution.