@Markiyan Hirnyk I gave you Brannan and Boyce as a general reference since it wasn't clear to me how much you knew about these things. But you are right, B&B doesn't cover the present situation.

You may want to take a look at

http://www.scholarpedia.org/article/Chetaev_function

After changing variables the present system can be written in vector form as

dx/dt = Dx+|x|e(x), where e(x) ->0 as x->0 and D is a diagonal matrix having the eigenvalues in the diagonal, say in the order -1, 0, 6.

As a Chetaev function you can take V(x1,x2,x3) = -x1^2 + 6*x3^2.

I leave it as an exercise for you to work out the details.

@Markiyan Hirnyk There are quite a few very good textbooks on the subject. The last one I used in an undergraduate class I taught was by Boyce and Brannan.

@Markiyan Hirnyk There are quite a few very good textbooks on the subject. The last one I used in an undergraduate class I taught was by Boyce and Brannan.

I see the point of using ApproximateInt since the integrals are replaced by finite sums and since W(r) and M(r) are just constants in the integration. I give a simple example below.

I notice, that vrad(x) depends on r as well as x, is that intended? Or should (1+r)^(-1/2) have been integrated over too? In your procedure

vproc1 := proc (X, Theta0, Beta, ThetaR) Wsol(parameters = ['theta0' = Theta0, 'beta' = Beta, 'thetaR' = ThetaR]); evalf[15](eval(vradsol, [x = X, 'theta0' = Theta0, 'beta' = Beta, 'thetaR' = ThetaR])) end proc

r doesn't appear.

The use of ApproximateInt seems OK to me. Here is an example.

restart;
eqI := diff(M(r), r) = Int(sqrt(x+M(r)), x = 0 .. M(r));
dsolve({eqI,M(0)=1},numeric); #Doesn't work
diff(eqI,r); #Doesn't seem possible to create a system not involving integrals
#Now ApproximateInt:
eqA := diff(M(r), r) =value( Student:-Calculus1:-ApproximateInt(sqrt(x+M(r)), x = 0 .. M(r),output=sum) );
dsolA:=dsolve({eqA,M(0)=1},numeric,output=listprocedure);
Mn:=subs(dsolA,M(r)):
RI:=subs(M(r)=Mn(r),rhs(eqI));
RA:=subs(M(r)=Mn(r),rhs(eqA)):
plot(RI-RA,r=0..1); #reasonably small

#################

However, to avoid using ApproximateInt you could (certainly in this example) use the procedural form of dsolve/numeric:

solproc := proc(N, r, Y, YP)
    YP[1] := evalf(Int(sqrt(x+Y[1]),x=0..Y[1]));
end proc;
dsol := dsolve(numeric, number=1, procedure=solproc,start=0, initial=Array([1]), procvars=[M(r)]);
plots:-odeplot(dsol,[r,M(r)],0..1);

I see the point of using ApproximateInt since the integrals are replaced by finite sums and since W(r) and M(r) are just constants in the integration. I give a simple example below.

I notice, that vrad(x) depends on r as well as x, is that intended? Or should (1+r)^(-1/2) have been integrated over too? In your procedure

vproc1 := proc (X, Theta0, Beta, ThetaR) Wsol(parameters = ['theta0' = Theta0, 'beta' = Beta, 'thetaR' = ThetaR]); evalf[15](eval(vradsol, [x = X, 'theta0' = Theta0, 'beta' = Beta, 'thetaR' = ThetaR])) end proc

r doesn't appear.

The use of ApproximateInt seems OK to me. Here is an example.

restart;
eqI := diff(M(r), r) = Int(sqrt(x+M(r)), x = 0 .. M(r));
dsolve({eqI,M(0)=1},numeric); #Doesn't work
diff(eqI,r); #Doesn't seem possible to create a system not involving integrals
#Now ApproximateInt:
eqA := diff(M(r), r) =value( Student:-Calculus1:-ApproximateInt(sqrt(x+M(r)), x = 0 .. M(r),output=sum) );
dsolA:=dsolve({eqA,M(0)=1},numeric,output=listprocedure);
Mn:=subs(dsolA,M(r)):
RI:=subs(M(r)=Mn(r),rhs(eqI));
RA:=subs(M(r)=Mn(r),rhs(eqA)):
plot(RI-RA,r=0..1); #reasonably small

#################

However, to avoid using ApproximateInt you could (certainly in this example) use the procedural form of dsolve/numeric:

solproc := proc(N, r, Y, YP)
    YP[1] := evalf(Int(sqrt(x+Y[1]),x=0..Y[1]));
end proc;
dsol := dsolve(numeric, number=1, procedure=solproc,start=0, initial=Array([1]), procvars=[M(r)]);
plots:-odeplot(dsol,[r,M(r)],0..1);

My value for k was chosen so that my dt agreed with yours. Their difference seemed to be zero judging from the plot, simplify didn't do it, but I was satisfied. However, as you will notice, I kept using yours.

I had fun making an animation of the jump and catch:

My value for k was chosen so that my dt agreed with yours. Their difference seemed to be zero judging from the plot, simplify didn't do it, but I was satisfied. However, as you will notice, I kept using yours.

I had fun making an animation of the jump and catch:

@brian bovril Thanks for the reference. I have uploaded a worksheet. It doesn't use the parameters option in dsolve/numeric, since I seem to be getting inconsistent results using that option in this problem. At this point I don't understand why results are inconsistent.

'implicitplot' takes much longer, but plotting is just used for guidance. 

MonkeyAndCoconut.mw

Addendum:

The inconsistency results in the parameters version apparently is due to fsolve remembering results.

Below the problem is shown (in p1) together with two solutions to the problem (p2 and p3).

But first the simple problem is solved analytically.

In p2 X is a local variable, in p3 X is global, but the statement forget(fsolve); is inserted before fsolve gets to work.

The parameters version will work fine also for the the original problem e.g. by inserting the statement forget(fsolve); as below in p3.

restart;
eq := diff(x(t), t) = a*x(t);
dsolve({eq, x(0) = 1});
solve(rhs(%)=2,t);
#############################################
restart;
eq := diff(x(t), t) = a*x(t);
dsol := dsolve({eq, x(0) = 1}, numeric, output = listprocedure,parameters=[a]);
X:=subs(dsol,x(t));
p1:=proc(a) local tp;
   if not type(a,realcons) then return 'procname(_passed)' end if;
   dsol(parameters=[a]);
   tp:=fsolve(X(t) = 2, t);
   if not type(tp,numeric) then undefined else tp end if
 end proc:
plot(p1(a),a=0..50);
p1(10);
p1(200);
#############################################
restart;
eq := diff(x(t), t) = a*x(t);
dsol := dsolve({eq, x(0) = 1}, numeric, output = listprocedure,parameters=[a]);
p2:=proc(a) local tp,X;
   if not type(a,realcons) then return 'procname(_passed)' end if;
   dsol(parameters=[a]);
   X:=subs(dsol,x(t));
   tp:=fsolve(X(t) = 2, t);
   if not type(tp,numeric) then undefined else tp end if
 end proc:
plot(p2(a),a=0..50);
p2(10);
p2(200);
restart;
eq := diff(x(t), t) = a*x(t);
dsol := dsolve({eq, x(0) = 1}, numeric, output = listprocedure,parameters=[a]);
X:=subs(dsol,x(t));
p3:=proc(a) local tp;
   if not type(a,realcons) then return 'procname(_passed)' end if;
   dsol(parameters=[a]);
   forget(fsolve);
   tp:=fsolve(X(t) = 2, t);
   if not type(tp,numeric) then undefined else tp end if
 end proc:
plot(p3(a),a=0..50);
p3(10);
p3(200);

@brian bovril Thanks for the reference. I have uploaded a worksheet. It doesn't use the parameters option in dsolve/numeric, since I seem to be getting inconsistent results using that option in this problem. At this point I don't understand why results are inconsistent.

'implicitplot' takes much longer, but plotting is just used for guidance. 

MonkeyAndCoconut.mw

Addendum:

The inconsistency results in the parameters version apparently is due to fsolve remembering results.

Below the problem is shown (in p1) together with two solutions to the problem (p2 and p3).

But first the simple problem is solved analytically.

In p2 X is a local variable, in p3 X is global, but the statement forget(fsolve); is inserted before fsolve gets to work.

The parameters version will work fine also for the the original problem e.g. by inserting the statement forget(fsolve); as below in p3.

restart;
eq := diff(x(t), t) = a*x(t);
dsolve({eq, x(0) = 1});
solve(rhs(%)=2,t);
#############################################
restart;
eq := diff(x(t), t) = a*x(t);
dsol := dsolve({eq, x(0) = 1}, numeric, output = listprocedure,parameters=[a]);
X:=subs(dsol,x(t));
p1:=proc(a) local tp;
   if not type(a,realcons) then return 'procname(_passed)' end if;
   dsol(parameters=[a]);
   tp:=fsolve(X(t) = 2, t);
   if not type(tp,numeric) then undefined else tp end if
 end proc:
plot(p1(a),a=0..50);
p1(10);
p1(200);
#############################################
restart;
eq := diff(x(t), t) = a*x(t);
dsol := dsolve({eq, x(0) = 1}, numeric, output = listprocedure,parameters=[a]);
p2:=proc(a) local tp,X;
   if not type(a,realcons) then return 'procname(_passed)' end if;
   dsol(parameters=[a]);
   X:=subs(dsol,x(t));
   tp:=fsolve(X(t) = 2, t);
   if not type(tp,numeric) then undefined else tp end if
 end proc:
plot(p2(a),a=0..50);
p2(10);
p2(200);
restart;
eq := diff(x(t), t) = a*x(t);
dsol := dsolve({eq, x(0) = 1}, numeric, output = listprocedure,parameters=[a]);
X:=subs(dsol,x(t));
p3:=proc(a) local tp;
   if not type(a,realcons) then return 'procname(_passed)' end if;
   dsol(parameters=[a]);
   forget(fsolve);
   tp:=fsolve(X(t) = 2, t);
   if not type(tp,numeric) then undefined else tp end if
 end proc:
plot(p3(a),a=0..50);
p3(10);
p3(200);

@brian bovril 

I don't know where the expression dt came from, nevertheless I tried the following.

restart;
eq1 := diff(y(t), t, t) = -49/5-0.1350000000e-2*Pi*sqrt((diff(x(t), t))^2+(diff(y(t), t))^2)*diff(y(t), t);
eq2 := diff(x(t), t, t) = -0.1350000000e-2*Pi*sqrt((diff(x(t), t))^2+(diff(y(t), t))^2)*diff(x(t), t);
dsol := dsolve({eq1, eq2, x(0) = 0, y(0) = 0, D(x)(0) = V*cos(theta),D(y)(0) = V*sin(theta)}, numeric, output = listprocedure,parameters=[V,theta]);
#The version above didn't work in Maple 12, which I used for my original answer, but it works in Maple 16.
X,Y:=op(subs(dsol,[x(t),y(t)]));
dt := -(25/126)*ln(-(-2+exp(-(324/1625)*Pi+(27/1625)*Pi*yp)+2*sqrt(1-exp(-(324/1625)*Pi+(27/1625)*Pi*yp)))/exp(-(324/1625)*Pi+(27/1625)*Pi*yp))*sqrt(78)/sqrt(Pi);
x0 := 45:
EQ:=proc(V,theta) local tp,drpt;
   if not type([V,theta],list(realcons)) then return 'procname(_passed)' end if;
   dsol(parameters=[V,theta]);
   tp:=fsolve(X(t) = x0, t);
   if not type(tp,numeric) then return undefined end if;
   drpt := evalf(eval(dt,yp = Y(tp)));
   tp-drpt
end proc:
EQ(25,1.3);
plots:-implicitplot(EQ(V,theta),V=0..50,theta=0..Pi/2);
th:=fsolve(EQ(40,theta)=0,theta=0.2..0.3);
dsol(parameters=[40,th]):
plots:-odeplot(dsol,[x(t),y(t)],t=0..2);

@brian bovril 

I don't know where the expression dt came from, nevertheless I tried the following.

restart;
eq1 := diff(y(t), t, t) = -49/5-0.1350000000e-2*Pi*sqrt((diff(x(t), t))^2+(diff(y(t), t))^2)*diff(y(t), t);
eq2 := diff(x(t), t, t) = -0.1350000000e-2*Pi*sqrt((diff(x(t), t))^2+(diff(y(t), t))^2)*diff(x(t), t);
dsol := dsolve({eq1, eq2, x(0) = 0, y(0) = 0, D(x)(0) = V*cos(theta),D(y)(0) = V*sin(theta)}, numeric, output = listprocedure,parameters=[V,theta]);
#The version above didn't work in Maple 12, which I used for my original answer, but it works in Maple 16.
X,Y:=op(subs(dsol,[x(t),y(t)]));
dt := -(25/126)*ln(-(-2+exp(-(324/1625)*Pi+(27/1625)*Pi*yp)+2*sqrt(1-exp(-(324/1625)*Pi+(27/1625)*Pi*yp)))/exp(-(324/1625)*Pi+(27/1625)*Pi*yp))*sqrt(78)/sqrt(Pi);
x0 := 45:
EQ:=proc(V,theta) local tp,drpt;
   if not type([V,theta],list(realcons)) then return 'procname(_passed)' end if;
   dsol(parameters=[V,theta]);
   tp:=fsolve(X(t) = x0, t);
   if not type(tp,numeric) then return undefined end if;
   drpt := evalf(eval(dt,yp = Y(tp)));
   tp-drpt
end proc:
EQ(25,1.3);
plots:-implicitplot(EQ(V,theta),V=0..50,theta=0..Pi/2);
th:=fsolve(EQ(40,theta)=0,theta=0.2..0.3);
dsol(parameters=[40,th]):
plots:-odeplot(dsol,[x(t),y(t)],t=0..2);

@Markiyan Hirnyk Please see the proof added above.

@Markiyan Hirnyk Please see the proof added above.

@J4James If alpha > 0 and beta > 0  then the first root found by Maple is real and the next two are imaginary. This follows from the fact that the cube root appearing in all the 3 results is the cube root of a positive number, thus real.

alpha*x+beta*x^3=f2;
r:=solve(%,x);
r[1];
op([1,1],indets(r[1],`^`(anything,1/3)));
#This is greater than or equal to
eval(%,alpha=0);
simplify(%) assuming real;
#which is >= 0.

#Of course if beta=0 and alpha > 0 there is only one solution.

If beta <0 you will observe that the situation is more complicated:

animate(complexplot,[[eval(r,{alpha=1})],f2=-10..10,thickness=3],beta=-5..5);

@J4James If alpha > 0 and beta > 0  then the first root found by Maple is real and the next two are imaginary. This follows from the fact that the cube root appearing in all the 3 results is the cube root of a positive number, thus real.

alpha*x+beta*x^3=f2;
r:=solve(%,x);
r[1];
op([1,1],indets(r[1],`^`(anything,1/3)));
#This is greater than or equal to
eval(%,alpha=0);
simplify(%) assuming real;
#which is >= 0.

#Of course if beta=0 and alpha > 0 there is only one solution.

If beta <0 you will observe that the situation is more complicated:

animate(complexplot,[[eval(r,{alpha=1})],f2=-10..10,thickness=3],beta=-5..5);