I have corrected the ODE.

You can select either the original definition of f() or f()->1 - just comment out the one you don't want

Rather obviously with f()->1, there will be no dependence on the variable 'R' in the final calculation

odeProb2.mw

Because when I run your worksheet I get no error from the PDEtools:-dchange() command.

One thing that looks odd is that that the picture of your code says PDEtools-dchange(), but the attached file says PDEtools:-dchange() - note the lack/existence of the semicolon

A further problem you are likely to have is that you define a quantity 'UPP' and then subsequently simplify a quantity 'UUP' which doesn't exist.

So on Maple 2017, the attached runs with no errors

dchangeProb.mw

If it doesn't run for you then specify Maple version and OS

There are so many errors in your latest worksheet that I have had to make many guesses about your intent.

As examples, I assume teat when your write sintheta, you actually mean sin(theta): and the same for cos(theta)

The attached is my best guess at what you intend, but I could be wrong! It calculates the integral you require for various values of the parameter 'R'. I have been unable to obtain aa analytic solution of your ODE, so  have had to resort to a numeric solution withg numeric values for 'R'

The final integral does not seem to vary very much as a function of 'R', wihc is ether true or I have committed a typo somewhere - if that is tru then ple\ase refence it in the attached, rather than producing another totally usless worksheet

odeProb.mw

Just for future reference - Maple (by default)  uses 200 points for plotting: so whatever the range you set along the x-axis, you start with 200points.

Now the the Mple plotting routines are "adaptive" so if they detect that the plot values are changing 'a lot' between successive plot points then they will insert more points. So the default of 200 is pretty much a minimum.

The  in-built adaptive point spacing *usually* works - but if you have a really 'spiky' plot (like yours) then sometimes the adaptive  routines fail, which is what happened in your case.

200 points (+adaptive) when the plot range is 100-1000 works.

200 points (+adaptive) when the plot range is 100-5000 doesn't.

It isn't so clever just to change the latter to 1000 points (+adaptive) - which is what setting numpoints=1000 does (more or less)

because there is a presupposition that user can actually read the relevant help, quoted below

Using Assumptions
In most cases, solve ignores assumptions on the variables for which it is solving.  However, the useassumptions option may be supplied to force solve to inspect any assumptions on the variables, and extract inequalities that it will add to the input.  For example:
solve(x^2-1,{x}, useassumptions) assuming x>0;
                            {x = 1}
is equivalent to:
solve({x^2-1, x>0},{x});
                            {x = 1}
but, in some case, the user may find the former more convenient The useassumptions option will not make use of any assumptions that cannot be rewritten as polynomial inequations.  In particular, it does not process assumptions restricting the domain:

Which part of this don't you understand??

I see no alternative.

I will submit an SCR, which may mean this bug will be fixed at some unspefied date in the futurw

If I use the Maple code (change the currentdir argument to somethng appropriate for your machine)

restart;
with(LinearAlgebra):
r_mean_Q_ini:=RandomMatrix(4,4):
r_mean_Q_ini_A:=RandomMatrix(4,4):
r_mean_Q_ini_B:=RandomMatrix(4,4):
legend := ["ini_que", "ini_que_A", "ini_que_B"]:
currentdir("C:/Users/TomLeslie/Desktop"):
ExportMatrix("test.mat",
             [r_mean_Q_ini, r_mean_Q_ini_A, r_mean_Q_ini_B],
             arraynames = legend,
             target = MATLAB
            );

restart;
currentdir("C:/Users/TomLeslie/Desktop"):
ImportMatrix("test.mat");

Then the ImportMatrix command, reimports the matrices with all names truncated to eight characters. However since the command returns a sequence of 2-element lists, with each 2-element list being [name matrix], at least I get all three matrices, even if two of them are associated with the same name.

However if I read this file into Matlab, then I only get two matrices, the first and the third. I'm assuming that since the second and third matrices are associated with the same name, then in some sense reading the third entry "overwrites" the second one.

I have played around with a variety of names for the matrices and the only conclusion I have come to, is that such names are always truncated to eight characters. Matlab accepts variable names with up to 63 characters. (Many, many years ago, I think this Matlab limit was 31.)

Hence I can think of no reason why Maple's ExportMatrix() would truncate names to eight characters - but it does:-(

No one here is going to retype your code.

Use the big green up-arrow in the Mapleprimes toolbar to upload a worksheet which demonstrates your problem

Question doesn't make much/any sense

Having obtained the numerical solution of the ODE for (probably incorrect) parameter values, it is relatively simple to use the output of this solution procedure to calculete any necessary function of (t, r(t), diff(r(t),t))

In the attached, I have shown one way to do this by calculating the integral given as eq(3.2) in your referenced paper.

odeProb2.mw

This seems to cover your requirement

"well how to find r(t) from your plot and integrate (3.2)"

However I have no idea what you expect me to do with

" we have a differential equation (3.3) ( i replaced theta by t ) , you solved numerically ( with your values for paramters M,l,C,Q,a,w) and find r against t diagram . now auther insert that in (3.2) then integrated  ( with theta_0=0.16 and r(0.159) ) and studied the relation between (3.2) and (2.6 T not t) that is figure 3 . "

because it doesn't make any sense to me. If you can define clearly what you want to calculate, I'm pretty sure I can do it

Consider

P:=(a,b,c,d,e,f)->a^2+b^2+c^2+d^2+e^2+f^2;

b:=x->x^2;
c:=x->x^2;
d:=x->x^2;
e:=x->x^2;
P(a,b(a),c(a),d(f),e(f),f)=0;

Now when you come to evaluate the last command in this group, the first thing that happens is that the arguments to P will be evaluated, so a-=a, b(a)=a^2, c(a)=a^2, d(f)=f^2, e(f)=f^2, f=f, so the last commant becomes

          P(a, a^2, a^2, f^2, f^2, f)=0,

so when this function is evaluated using the definition of P() above, this will become

a^2+a^4+a^4+f^4+f^4+f^2=0

and Maple correctly return a solution for 'a' in terms of 'f'

On the other hand, if you solve

P(a,b,c,d,e,f)=0,

then from the definition of P(), you will be solving

a^2+b^2+c^2+d^2+e^2+f^2=0

Once again Maple will return a solution for the equation you supply, in terms of the variables you specify

because I can get a usable definition for 'f(r)', and some of the parameters. However I can't easily figure out values for the others, and I have no plans to read the whole paper in order to determine what these might be - so I just guessed a few values.

Maple can now produce a solution but it doesn't look much like the curves in the paper - so I assume my guesses on the 'unknown' parameters are way off. See the attached - and read the comments

odeProb.mw

because your ODE contains two dependent functions, r(t) and f(r), (which maybe should be f( r(t) )? ). Two unknown dependent functions require two ODEs, along with with a sufficient number of boundary/initial conditions (which will be at least one more than you have now).

Although it is sometimes possible to numerically solve IVPs in containing parameters, it is generally a better idea to have numerical values for everything, so a value for 'c' in the boundary condition would also be a good idea

As I stated in my earlier response (emphasis added)

I can however envisage situations where thsi symbolic approach might/would fail - probably not applicable to your "relatively simple" objective function. But consider, for example, the case of a non-continuous (hence non-differentiable) objective function - or possibly a periodic objective function where simple differentiation might produce an infinte number of possible extrema...

So just for amusement, I changed the definition of 'f' in your code to sin(x)*sin(y) - ie a 'somewhat' periodic function. Guess what - your code *fails*. 'Failure' here means that  I terminated it after about 5 minutes because it didn't seem to be getting anywhere. Don't feel too bad about this! I tried the same objective function with vv's code, with the same outcome.

So both you (and vv) have produced code which 'works' for a very, very, restricted class of objective functions. But as a general method of producing 'exrema' for *all* functions in a restricted domain?? Get Real!

Just for comparison, when I use the DirectSearch package on the same objective function, ie sin(x)*sin(y), with the same constraints, after about 10secs I get the output

Matrix(9, 3, [[1., [x = -4.71238898030787, y = -4.71238898010830], 77],
                    [1., [x = -1.57079632760966, y = -1.57079632159199], 77],
                    [1., [x = 1.57079632678335, y = -4.71238898042640], 53],
                    [.989983390184016, [x = -1.44523767572349, y = 4.64662778819133], 130],
                    [.948550029128931, [x = 1.29002111738053, y = 1.41067657599589], 145],
                    [.582556164299749, [x = 4.81840612640643, y = -6.90910816356876], 124],
                    [.229480375757100, [x = -3.52568245328446, y = .659108631202517], 127],
                    [.199812680179242, [x = 3.50512265552320, y = -2.54486188486310], 104],
                    [.151561851824201, [x = -1.14284786480444, y = -6.45055016424005], 78]])

ie three 'global' maxima within the specified domain, along with six 'local' maxima.

So serious question which is the *better*  code to use for a solution? Yours?, VV's?, DirectSearch()?

You state

I won't be able to run DirectSearch instructions if the packet isn't installed in the computers at college.

Have a friendly chat with those who do system administration at your college. This package is completely free: Download/installation would take less than 5mins - why wouldn't this be done???

I've just tried this for a hemisphere.

The first method in the worksheet executed, but the geodesic produced seems to be incorrect

The second method using CalculusOfVariations/EulerLagrange still works

See attached

geodFix2.mw

If I were going to study geodesics seriously, then I would probably rewrite the second method in this problem in both cylindrical and spherical polars. The former would handle geodesics on cones and cylinders; the latter would handle geodesics on spheres. I'm pretty sure that using the 'appropriate' coordinate system would make the ambiguity problem go away