@acer That's an incredibly complex code which does incredibly complex things!  It does the expected:

> 2*sin(x)*cos(x);
                          2 sin(x) cos(x)
> H(%);
                               sin(2 x)

and also the unexpected:

 > 4*sin(x)^2*cos(x);
                                       2
                               4 sin(x)  cos(x)

> H(%);
                               cos(x) - cos(3 x)

I would have expected 2*sin(2*x)*sin(x) as the answer but the combine step in the definition of H intervenes and, for better or worse, takes it one step further.

@Carl Love Your combine/half_angle code is very clever and it rightly belongs to Maple's default collection of convert(...) functions.  Until then, I will keep it as an add-on in my initialization file.

As to the double_angle conversion, applying combine is not ideal. For instance

expr := 4*sin(x)*cos(x)*sin(y)*cos(y);
                     expr := 4 sin(x) cos(x) sin(y) cos(y)
combine(expr);
                   1/2 cos(-2 y + 2 x) - 1/2 cos(2 y + 2 x)

I would have preferred the answer of sin(2*x)*sin(2*y). I don't know how to construct an extension to convert that will do that.

@Carl Love Thank you very much for that example. Using subsindets is new to me and I found your example very instructive.

Here is a generalized version of your construction. I intend to place in my Maple initialization file.

Download convert-to-half-angle.mw

I also played around with the idea of making something similar for the double-angle identities:

sin(x)*cos(x) = sin(2*x) / 2;
cos(x)^2 = (1 + cos(2*x)) / 2;
sin(x)^2 = (1 - cos(2*x)) / 2;

but did not get very far.  If you have suggestions about how to go about this, I will be eager to hear them.

@acer Thanks for the explanation.  I see the fault of my original approach and I sorry for having called it a "bug".

@acer Thanks for your detailed explanation. It's clear to me now why my approach wasn't working. It's the map[2] command that I needed. Thanks again.

@Kitonum The others have pointed out what I was doing wrong and how to fix it.  That said, I very much like your approach because that we we may do:

subsindets(sin(a+b), sin(anything),t->2*sin(op(1,t)/2)*cos(op(1,t)/2));

and obtain 2 sin(a/2 + b/2) cos(a/2 + b/2)..

That brings up a question:  Is is possible to apply two transformation rules in a single call to subsindets?  Specifically, can we call subsindents just once in order to apply the two rules defined in half_angle_rule in my original post to the expression sin(x)+cos(x)?

@Carl Love Thanks for observation. Now I see what I was doing wrong.

Both examples run fine on Linux.  I tried them in GUI and also in a character-based terminal.

Since Mac is really a unix clone under the hood, your error may stem from the stacksize limit in your unix shell.  I have no experience with Macs, so I don't know what sort of a shell is its defaul.  If it is something like the classical csh or tcsh, type "limit" at a terminal to see what your stacksize is. Mine says "stacksize 8192 kbytes".  You may raise the stacksize if needed.  If your shell is something like sh or bash, the command is "ulimit".

@ogunmiloro You have a linear system of 11 differential equations in 11 unknowns. For the equilibrium to be stable, all eigenvalues of the coefficient matrix should be on the left-hand side of the complex plane. The Routh-Hurwitz criterion (look it up in Wikipedia) provides necessary and sufficient conditions for that. 

In the attached worksheet I have extended my previous calculations to analyze the locations of the eigenvalues through the Routh-Hurwitz criterion. In the end, we see that the stability of the system is determined by two easy-to-verify explicit inequalities that should hold among the system's parameters. We also see that the inequalities may or may not hold, depending of the relative magnitudes of those  parameters.

Have a look at those inequalities and see whether you have enough information about your model to justify their validity.  In the worst case, you will need to supply numerical values for the parameters to verify that the inequalities hold.

This site refuses me to display the contents of my worksheet.  You will need to download it and view it in Maple.

Download worksheet: mw2.mw

The entries in your K matrix are large numbers such as 1.6123642440×10¹⁶. To handle such large numbers you are forced to increase Digits. As a result, Maple has to switch from hardware to software floating point calculations which greatly slows down your work.

Is there a way around this?  Sure!  Consider that you may express the distance from the Earth to the Sun as 150×10⁶ kilometres, or 150×10¹² millimetres, or 8 light-minutes, or 1 astronomical unit.  They all mean the same thing, but calculating with the distance set to 1 results in a much more stable numerical calculations than setting it to 150×10¹².

Consider changing the units of measurements in your model so that the coefficients reduce to a reasonable range, ideally things that are comparable to 1. Then you won't have to increase Digits and all the calculations will be done in hardware, speeding things up tremendously. There is nothing sacred about the use of meter or kilogram for the units of measure. Pick you own units.

You are likely to receive better responses if you provide a specific example that you are interested in. Examples made up by others may not do what you want.

The problem you have outlined is obviously too complex for a first attempt.  I suggest that you begin with a simpler problem in order to get used to both the mathematics of the problem, and the corresponding implementation in Maple.

Begin with finding out about (1) forward Euler scheme versus backward Euler scheme; (2) The CFL (Courant-Friedrich-Lewy) stability condition.  [These are mathematical concepts, not Maple.]

Once you have understood the above concepts, try solving the simplest diffusion equation:

There is no point in struggling with your original problem if you cannot solve this significantly easier one.

@JAMET Near the top of the code, change L := 5*R to L := 3*R.  And also change a=1.3*H to 1.8*H.

If that doesn't do what you want, and you want further input, then please post your maple worksheet (the *.mw file) as I had asked before.  I won't be able to spend any more time with the kind of code that you are posting.

@JAMET You may add any object that you wish to the skeleton.  Here I have added a rather arbitrary polygon,  Edit/change as needed.

In the worksheet I have highlighted the changes (which are very few) relative to the previous "skeleton" version.

Worksheet: pump-head.mw