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.

I have fixed several issues with you worksheet as noted.  See if this does what you want.

Download mw.mw

Maple's Explore is the right tool for solving your problem.  Explore assigns graphical "sliders" to the problem's parameters.  You may change the parameter values through the sliders and observe the effect on the problem's solutions.

Your worksheet contains seven equations, while you said you are interested in solving three of them,  I have extracted those three equations into a new worksheet and have added the Explore command.  I have assumed that the parameters take values in the range -5 to +5.  You may change those ranges as needed.

Here is the worksheet:  explore.mw

In your worksheet you are attempting to solve a system of six differential equations by splitting the process into two stages.  In the first stage you solve for the unknowns p, q, r, and then in the second stage you solve for theta, phi and psi which depend on p, q, r.

The preferred way is to solve the entire system of six equations in one pass.  I have modified your worksheet to do that.  Along the way I have made several other minor adjustments.  I believe these should be acceptable to your Maple 13, Let us know if it doesn't work.

Пример02-новый.mw