@mmcdara 

Thank you so much for your time and intervention. Your comments were particularly useful, particularly "In fact the condition Norm(err) <= tol is never met." I applied your comments and adjusted the tolerance, the codes now works perfectly.

Once again, Thank you.

@mmcdara 

Thank you so much for your observation. it is appreciated. I totally agreed with your submission that h is nor increasing. The motivation for writing 

h := 0.1 * h * (c + 1) * (tol/Norm(err))^(0.2)

is to ensure that no fixed step exists during the computation rather a variable one. To ensure this, a formular given by

h[new]=k*h[old]*(tol/Norm(err))^(1/(order+2)) in this case, order=3 is used.

Honestly, your time and comments are appreciated.

@mmcdara 

Thank you so much for your time and the modification made. It is appreciated.

However, I am yet to obtain the desired result. I observed from the results that the value of h is increasing geometrically instead of decreasing. I believe this may account for not obtaining the desired result.

Please, it would be appreciated if you could help in this regards.

Thank you, and best regards.

@acer I want a 3 dimension figure that connects the time (step sizes), the exact values and the errors column.

Thank you for your time.

@rcorless 

Thank you so much for your comments. You have indeed opened my eyes to see beyond and motivated me to want to learn more. It is indeed appreciated. Permit me to quote some of your statements then I reply

1. "Another way to modify your code is to write your own numerical boundary-value problem solver using (for instance) quasilinearization (your problem is a second-order nonlinear boundary-value problem, and the process of quasilinearization replaces nonlinear problems with a sequence of linear ones --- see eg https://en.wikipedia.org/wiki/Draft:Quasilinearization for a description).  Is that what you want to do?"

This is indeed what I want. Ability to write my own numerical boundary-value problem solver. Please can you be of help? Also, if possible can I learn from you?

2. Your code looks as if it might be a fixed stepsize finite-difference approach.

You are 100% right.

3. Are you trying to learn how "shooting" methods work?

Would like to know this also.

Thank you once again for your time and valuable comments.

@tomleslie Thank you for your comment and observation.

Actually, the correct values occupy the even positions in columns Ex.y, Ex.z, Error y, and Error z respectively. Also, could there be a reason why consecutive values in columns Ex.y and Ex.z always come in "pairs"? It would be appreciated if the code can be modified so that the pairing values will not occur.

Thank you for your time and best regards.

@acer Thank you so much for the enlightenment. It is appreciated indeed and thank you for your time

@mmcdara Thank you for your time and response. You right. Please kindly check the attached document above for your perusal. Your comment would highly be appreciated.

@tomleslie 

Thank you for your time and response. However, there exist differences. See the attached. I really wish to know what is responsible for the differences and when does one use the syntaxes. Also, which is the best of the two

SDFFM_STAB.mw

@dharr

Thank you so much. It is appreciated 

@dharr 

Thank you so much for your kind reply. Is it possible to remodify the code? If yes, can you be of help?

Thank you so much and kind regards.

@Preben Alsholm Thank you so much. I didn't take note of that. Your comment is appreciated. 

@Carl Love Your guess is correct. the two numbers were fictitious numbers to complete the list in the vector.

@tomleslie I really appreciate your suggestion and time. However, when I run the code in Maple 2016, there is an error message "Error, (in plots:-display) display does not accept the legend option". I have tried to check what may be wrong but no luck. Could you please be of help?

@Carl Love Thank you so much. Your comments and suggestions were useful. It is appreciated