@Markiyan Hirnyk It appears that VectorCalculus:-int with the Rectangle input returns an iterated integral.
Let phi be unassigned:

VectorCalculus:-int(phi(x,y), [y, x] = Rectangle(-infinity .. infinity, -infinity .. infinity));
lprint(%);
int(int(phi(x, y), y = -infinity .. infinity), x = -infinity .. infinity)

I think this explains why VectorCalculus:-int with Rectangle works correctly for the function g. In my answer I used iterated integrals instead of double range integrals to get the correct answer.

I simply don't understand what you are talking about. Could you be more specific? Has it got to do with Maple?

@Christopher2222 Quite right; you can't.

@shakuntala In the sense I described above the ode involving only f (and with the boundary conditions with theta left out) does indeed have two solutions (at least). I only looked at the first value in B, i.e. B[1] = 1.
The other solution for f is negative and decreasing. I'm not familiar with the physics of this, but this may be an unwanted solution for those reasons alone.
If not, then there is a mathematical (or at least a numerical) problem for the total system with that solution: The numerator in the theta equation when solved for diff(theta(eta),eta) is 7*f(eta)+30:

That is bad news because 7*f(eta)+30 becomes zero in the interval 0..1.1 as seen in this graph:

@shakuntala Do you have just one solution? Why talk about dual solutions if you don't even have one?

## But to clear up matters: Do we mean the same by the expression dual solution?

In my use of the term "dual solution" it means that a given bvp problem has two solutions, as e.g. the simple algebraic equation x^2 = 4 has two solutions (-2 and 2).

@shakuntala Plotting two (or more) solutions is no problem at all. In fact it seems that your worksheet has the code for it.
Your problem (as it has been presented so far) is to find any solution at all.

Have you done that?

You wanted to impose two conditions on theta and three on f, yet the total order of the system is 3+1=4, as I pointed out. You must leave out one of your five conditions!
You can e.g. leave out theta(1.1) = 0 or you can leave out D(f)(1.1) = 0 (just examples, but in both cases dsolve gives results without problem.

@vv I tried an animation of a disk with a "background" consisting of your 4 plots (p1,p2,p3,p4) and a texplot.
It just confirms what has already been described by you.
 

restart;
with(plots):with(plottools):
p1:=plot( x^3,x=-1..1,thickness=20,color=red):
p2:=plot(-x^3,x=-1..1,thickness=20,color=blue):
p3:=display(rectangle([0.5, 1],[0.75,-1],color=green)):
p4:=display(rectangle([-1,-0.1],[1,0.1],color=yellow)):
#######################################################
pt:=textplot([0,0.5,"Happy New Year", 'font'=["times","roman",20]]):
Q:=a->display(disk([0,0],a,color="NavajoWhite"),axes=none):
plots:-animate(Q,[a],a=0..1.25,background=[p1,p2,p3,p4,pt],paraminfo=false);

You don't explain what plrPt is.
In order to transform the 2D animation to a 3D animation you would have to transform each frame.

@Christopher2222 To combine your animated field with the fixed field you can produce both and display them together:

with(DETools):
LVS := [diff(x(t),t)=x(t)*(1-y(t)),diff(y(t),t)=.3*y(t)*(x(t)-1)];
DEplot( LVS, [x(t),y(t)], t=0..12, x=0..3, y=0..2, arrows=comet,
  animatefield=true,numframes=30,dirfield=50); p1:=%:
DEplot( LVS, [x(t),y(t)], t=0..12, x=0..3, y=0..2,color=gray); p2:=%:
plots:-display(p1,p2);

But here you have just one ode. In Kitonum's revised example there were 90 odes, one for each frame.

@Christopher2222 The field does indeed move in Kitonums animation (New version).

If you add animatefield you get an animation having each frame an animation. Just like having a movie whose individual frames are movies themselves. How do you expect that to play?
The movie could be it all of Hitchcock's movies played in sequence chronologically according to release date?

@Christopher2222 Adding animatefield=true to the DEplot command makes an animation for each of the 90 frames in the outer animation. I don't think that is what you intended.

@Christopher2222 DEplot accepts the option maxfun= ... as does dsolve/numeric. Try maxfun=10^6 or whatever you think will work. In which setup are you using animatefield, i.e. which code?

@asa12 animate in Maple versions 12, 13, and 14 (but no later ones) is using subs in places where it ought to use eval.
## I now see in a reply to Kitonum that you actually use Maple 15. There my code works.
## Be aware that if you have axes = none, then what you see initially is just white space.
## Since the following may be useful to someone using earlier versions, I shall leave it as is.
##
One solution is to redefine `plots/animate` as follows.

restart;
## The first line redefines animate:
`plots/animate`:=subs(subs=((x,y)->eval(y,x)),eval(`plots/animate`)):
## Now test it:
sols := [[0,0], [6,0], [6,4], [0,4], [3,6], [6,4], [0,0], [0,4], [6,0]];
plots:-animate(plot,[sols[1..round(i)],thickness=3],i=1..nops(sols),paraminfo=false);

I just tried in Maple 12, and it works for me. I used this device when Maple versions earlier than 15 were the current ones.
As an alternative you can use this two step method:

restart;
sols := [[0,0], [6,0], [6,4], [0,4], [3,6], [6,4], [0,0], [0,4], [6,0]]; 
p:=proc(i) plot(sols[1..round(i)],_rest) end proc;
## Test the procedure p:
p(4,color=blue,thickness=3);
## Now animate p:
plots:-animate(p,[i,thickness=3,axes=none],i=1..nops(sols),paraminfo=false);

Using display/insequence as Kitonum does is perfectly fine. My point was just that animate can be used, if you want to.

Here is a procedure that parametrizes a polygon defined by a list of points given as [x,y].
(Since first posting this I changed it slighly so it can be used in any dimension).

par:=proc(pts::listlist,t::name) local d,n,i,j,r,s,p,px,py;
  d:=nops(pts[1]);
  n:=nops(pts);
  for i from 1 to n do
    j:=modp(i,n);
    r[i]:=pts[i]*(1-s)+pts[j+1]*s
  end do;
  p:=piecewise(op(op~([seq([t<i,eval(r[i],s=t-(i-1))],i=1..n)])),undefined);
  seq(evalindets(p,list,L->L[i]),i=1..d)
end proc;
## Now using CS as defined above we can do: 
px,py:=par(CS(12),t): 
plots:-animate(plot,[[px,py,t=0..T],thickness=7,color="SkyBlue"],
T=0..nops(CS(12)),paraminfo=false,axes=none,frames=100);

Here is an example of a space curve:
 

r:=rand(1..5.0);
L:=[seq([r(),r(),r()],i=1..100)]:
px,py,pz:=par(L,t):
plots:-spacecurve([px,py,pz],t=0..100,thickness=5,numpoints=5000);
plots:-animate(plots:-spacecurve,[[px,py,pz],t=0..T,thickness=5,numpoints=5000],T=0..100,frames=500);

Since it is Christmas here is a version with stars:
 

restart;
CS:=n->[seq((3+(-1)^k)*[cos(Pi/n*k),sin(Pi/n*k)],k=0..2*n)];
N:=8:
plots:-display(Array([seq(plots:-animate(plot,[CS(n)[1..round(i)],thickness=3,axes=none],i=1..nops(CS(n)),paraminfo=false),n=3..N)]),scaling=constrained);

It seems that I cannot show a gif file of the whole array, but I can give you one of the animations.
Here is CS(5):
 

plots:-animate(plot,[CS(5)[1..round(i)],thickness=3,axes=none],i=1..nops(CS(5)),paraminfo=false);