@Earl I don't see an easy way to modify the code to do what you have asked.

There is some literature on the subject.  Perhaps a good way is to Google for "how to calculate the geodesic going through two given points".  I just did that and obtained two promising articles.

1. This one appears to have been written by undergraduates.  Their methods seem to work but are computationally inefficient.
2. This one appears to be more sophisticated. I haven't looked at the details but it seems to be worth a look.

I may give this question a try when I have some more time but not now.

@ecterrab That works!. Thanks!

PS: In the first of the two maplinit lines that you have provided, the "=" should be ":=".  

@ecterrab Thanks for the quick fix, as always!

Edgardo, that's a very nice enhancement.

Would it be possible to let users make the subscripted variables the default by setting some kind of environment variable in mapleinit?

@vv Yes, pdsolve() does return the correct general solution.  Something goes wrong when it tries to apply the initial condition.

Here is a possible enhancement to Maple pdsolve().  In the same way that dsolve() offers an implicit option. it would be helpful if pdsolve() also offered that option.  Then the solution of the initial value problem u_t + u u_x = 0, u(x,0) = f(x) would be expressed as f(x − ut) = u which is more pleasant than the RootOf expressions to look at. 

@Earl That looks good.  Good job!

@Preben Alsholm I like that.  Very nice!

@acer Thanks for getting to the bottom of this.  Now it's clear to me why D[1](U)(0,t) does not return a value .

So now we have two ways for calculating the boundary flux.  One way is through the boundary_flux() function which I described earlier.  Another way is to calculate −D[1](U,h,t) where h=Float(1,−Digits/2).  That actually calculates the flux at y=h, but since h is so small, the value calculated here is practically the same as the flux at the boundary.

As long as we are at it, let me point out that a more accurate version of the boundary_flux() function may be obtained by replacing the 1st order finite difference by a 2nd order finite difference, as in:

h := Float(1,-Digits/2);
boundary_flux := t -> - ( -3/2*U(0, t) + 2*U(h,t) - 1/2*U(2*h,t) ) / h;

This would be my preferred method of calculating the boundary flux if I had a need for it in my work.

@Moqifang I don't know what prevents the evaluation of the y derivative at the boundary.  It's not a version issue, since it doesn't work on my Maple 2022 either.  Nevertheless, there is an easy workaround as we see in the following fragment corresponding to your screenshot:

@wmcnally I see that in the case of the pendulum equations you have

sys := A:-Linearize(format=all)

while in the case of the RLC circuit you also specify params and paramcheck.  Perhaps that's the reason for the pendulum case not working?  That's just an observation; I know nothing about Maplesim.

If you don't receive further suggestions to your question here, you may want to ask for help by emailing <support@maplesoft.com>.

@Rakshak I don't know how to apply the reduction automatically to every expression.  Sorry.

I don't have Maplesim, so I cannot view your worksheet.  Did Maplesim calculate those linearized equations or did you?  In any case, those equations are not correct.  Here are the linearized equations of the double-pendulum as calculated in Maple:

or equivalently:

where I have written u(t) and v(t) for the angles to simplify the typing.

@tomleslie You are bringing up multiple and mutually disjoint issues.

1. Mathematical modeling issue.  You note that a more accurate model will account for a thermal impedance at a contact area between two bodies of differing temperatures.  Okay, one may take that into account, but then we would be solving a different problem, not the one that the OP has asked.
   
   Even if we account for the impedance, the model may still be criticized as "physically inaccurate" since it ignores that everything is made of atoms and the model is ignoring those. Where does one stop?
2. Theoretical considerations.  Within the context of the original model, we have well-established existence and uniqueness theorems of OP's classical initial/boundary value problem.  A unique solution exists even when boundary conditions are discontinuous. In fact, the discontinuities need not be only at the "corners".  A discontinuous initial condition

   u(y,0) = piecewise(y< N/2, 0, 1);   #  0 < y < N

   is perfectly acceptable as well.  In all cases we have an infinitely smooth solution away from the boundaries, but the solution may have discontinuities at the boundaries.  Nothing wrong with that.

3. Numerical solution. As noted in the previous item, the problem has a well-defined unique solution.  Whether a numerical algorithm can capture that unique solution is a different issue. As you have noted, tweaking the parameters yields various candidates for the solution, each of which is an attempt, by the designer of the numerical algorithm, to capture the problem's one and only true solution. The failure of a numerical algorithm should not be construed as a shortcoming of the mathematical model as you seem to be implying.

@tomleslie The mismatch between the boundary values at (0,0) is not a problem.  In fact, that's a very common occurrence and is harmless as far as solving the PDE is concerned.  All it says is that the solution is discontinuous at that point on the boundary.

Consider the following scenario which leads to such a discontinuity.  Take a metal bar which is initially at temperature u=0 throughout.  At time t=0 bring one end of the bar in contact with the wall of a hot (say 300 degrees) oven.  The temperature in the bar will evolve according to the heat equation u_t = u_xx.  The initial condition is u(x,0)=0.  The boundary condition at x=0 is u(0,t)=300.  The temperature has a discontinuity at (0,0) but the bar doesn't mind.

@acer Very clever idea.  Vote up!