I assume there are points on the border that are placed in the same way as the red dots, aren't there?

If it is so and you replace the origin by a quarter of circle with a small radius,  then there exists an homeomorphism from this truncated domain to a square (more simply you can unfold and rescale the truncated domain to get a square).
So you find yourself in the situation you addressed in the example contained in your worksheet.

In fact not really the same because the origin (more precisely circle C) has been replaced by a side of the square which thus contains N boundary points.
So, after having decomposed the matrix the way you did, simpliy remove on each matrices N-1 rows and columns corresponding to N-1 of these boundary points .

It seems to me that as long as:

1. the points are vertices of a structured mesh made of squares (possibly curved squares),
2. and there exists an homeoporphism from the domain you consider into a square (possibly after having removed some singularities),

then you can apply your decomposition strategy.

Note that if your domain was the whole disk D, replacing the origin by a small radius circle would create a fictitious boundary containing N points (there still exists an homeomorphism from D\{0} into the unit square), then you should remove N rows and columns in the matrices the decomposition gives.

@Ahmed111 

You're wrong, add these commands at the end of the file I sent you to see that "Maple DOESN'T treat both I and J as the imaginary unit"  (at least with my Maple 2015 version [see top this same file])

is(I=J)
                          false

sqrt(J);
sqrt(I);
                     1  (1/2)   1    (1/2)
                     - 2      + - I 2     
                     2          2         
                              (1/2)
                             I     
J*conjugate(J);
I*conjugate(I)
                               1
                                _
                              I I

Because I use Maple 2015 and the simplification rules will not do the same thing that those of your Maple version, I did not go further 
the expression of f2.
Nevertheless you will found how to "factor out I = sqrt(-1) from square roots".
up_to_f2.mw

Use this in your own code and Maple version to see if this helps.

@Alfred_F 

You write "the three circles aren't enough to define exactly one sphere. There are obviously a whole host/set of spheres that meet the requirements": so why did you ask the question?


"This leads to the question: What is the geometric locus of all sphere centers?"
Ill posed problem: radii of the intersection circles doesn't enable finding this locus.

To understand why ;simply take the anologuous 2D situation of a disk D (it will be simpler to consider a disk than a circle) intersecting axes Ox and Oy.
Its imprints are two segments Sx and Sy whose mid points have coordinates (a, 0) and (b, 0). Thus the center of D is (a, b).
Let Lx and Ly the lentghts of 2_*p and 2_*q respectively. Then the square radius of D is R² = b² + p² = a² + q² [eq1]
Can we move Sx along Ox and Sy along Sy to other locations, let's say (a_*m, 0) and (b_*n, 0) in order that they still are imprints of a disk D' of a possibly different radius R'?
Because of the relation b²n² + p² = a²m² + q² [eq2] not all couple (m, n) is allowable.
Using [eq1], [eq2] rewritties n² = 1+ (m² -1)_*(p/q)²  [eq3].

One special case if p = q which leads to m = n or  m = -n.
Then the locus of the disk centers are the two main diagonals y = x and y = -x.

Let us assume that p < q, a > p and b > q (the reasoning which comes can easily be extend to the case p > q by interverting p and q), just to locate the center of disk D in the first quadrant.
When m = 0 the segment Sx corresponds to the interval x = [-p, p] and by relation [eq3] the vertical imprint Sy of the disk is the seqment [Z-q, Z+q] where Z is the positive square root of 1- (p/q)² (which is strictly larger than 0). In this case the radius of the disk is 2_*q (Z-q < 0 and Z+q > 0).
So the point [0, Z > 0] is on the locus if the disks centers.
The simple change of coordinates x --> x, y --> (y-Z)_*(p/q) shows that the locus of the disk centers is made of two two main diagonals (a point established before).
So, in the original referential the locus is made of two straight lines passing through the point (0, Z) and of slopes arctan(b/a) and -arctan(b/a).

Note that the sole values of the imprint radii p and q are not enough to determine this locus.

Extending this reasonning to the the sphere case  shows the locii are straight lines (the 4 main diagonals of the 8 octants of the ad hoc transformed referential).
This transformation is of the form x --> (x-a_x)/p_x, y --> (y-a_y)/p_y, z --> (z-a_z)/p_z, where the a_i denote the coordinates of the imprints  mid point and the p_i the imprint sizes.
By the inverse transformation one gets the equations of the locci.

Lets take a sphere S of radius r whose center C has coordinates (a, a, a).
So the intersections of S with each of the three canonic planes are circles of identical  radii (provided of course that the intersections do exist, that is that r > a). Let u this common radius
Now move the center at  another place (A, A, A) (A > a), rgere is always a value or R > r such that R²-A²=u².
So the radius and the coordinates of the sphere cannot be uniqueky determined by the radii of the intersection circles.

"the sphere's radius r and the coordinates of its center can be determined?": NO

One has exactly the same result in 2D for a circle intersecting two orthogonal lines.

@KIRAN SAJJAN 

saying "When I change the value of j, only one line is plotted".
If you make a zoom you will see there are 4 curves. Change for instance Grt from 0.1 to 2 to convince yourself the code gives different resullts.
For instance

By the way I wonder by what mistery you can have past 12 hours working on the code I posted only 8 hours ago?

@KIRAN SAJJAN 

I only took Br=0.1 has an illustration of what you can do for the other values in BrVals.

Feel free to go further when and proceed as described in the second part of the worksheet for those other Br values.
I don't think it is an impossible task for you, you have just to change the value of j here
As I told you I don't have more time to spend on this problem and your initial worksheet symetry_paper_work.mw proves you know perfectly how to browse all Br values in BrVals, build a table of results and display them...so I don't feel you need me any longer.

@KIRAN SAJJAN 

I guess this is closer to what you have in mind:

• The first part of the worksheet 2_domains_sand15 explains the strategy I use.
  The whole solution (the "two-domains" solution) depends on two parameters A and B where A is the common unknown value of u₁(0) and u₂(0) and B is the common unknown value of Theta₁(0) and Theta₂(0).
  In domain D₁ = {y, y in [-1, 0]} the differential system is well posed and can be solved numerically a soon as numerical values a and b have been assigned to A and B. Let S₁(y ; a, b) the solution of the problem 1 on domain D₁.
  In domain D₂ = {y, y in [0, +1]} the differential system can also be solved numerically and let S₂(y ; a, b) its solution.
  Let Jump the vector whose components are

  Jump1 := (D(u1))(0)     - Z1*(D(u2))(0),
  Jump2 := (D(Theta1))(0) - Z4*(D(Theta2))(0)
  

  For an arbitrary choice of a and b there is no reason for Jump₁ = Jump₂ = 0  when S₁(0 ; a, b) and S₂(0 ; a, b) are used to assess D(u₁))(0), D(u₂))(0), D(Theta₁)(0) and D(Theta₂)(0).
   
• Second part of the worksheet 2_domains_sand15.
  I wrote a procedure named TwoDomainsSolver(a, b) wich gathers all the steps presented in the first part and returns the L₂ norm of the Jump vector.
  Then all you have to do is to use an optimization function to findthe couple (a*, b*) wich minimizes the output of TwoDomainsSolver(a, b) over some (a, b) domain.
  Once (a*, b*) is found you can draw the "left" and "right" solutions.

I provide you two worksheets corresponding to two different data
2_domains_sand15.mw
2_domains_sand15_other_data.mw

For other data, for instance higher values of the Grashof number, you may get this error

Error, (in dsolve/numeric/bvp) initial Newton iteration is not converging

A quite common one but an error which is always borring and may be relatively difficult to fix, even if some recipes are given here Numerical Solution of Difficult ODE Boundary Value Problems

One example is given in 2_domains_sand15_Gr=15.mw where this error appears for Gr=15.
You will find in this worksheet a way to fix it (read carefully the last comment lines).

I'm sory to tell you that I'm done with this thread because I am a lot of other things to do.
I hope this response has been helpful.

@KIRAN SAJJAN 

Your initial problem contains:

• 4 bc given at point y=-1 and/otr point y=+1
• 4 continuity equations written at point y=0

This cannot ne solve this way and a "true" nv problem requires 8 bc written at points y=-1 and y=+1.
So I rewrite your initial problem in term of a well-posed bvp but introducing some "degrees of freedom":

u2(-1)     = au, 
Theta2(-1) = aT, 
u1(1)      = bu, 
Theta1(1)  = bT,

where au, aT, bu and bT are formal quantities and not numbers.

Maple does not allow solving a bvp with "parameters" (as it does for initial value problems).
So I defined a procedure F(au, aT, bu, bT) which takes as input 4 numbers and deliver a measure of the error E on the continuity equations

u1(0) - u2(0) = 0,
Theta1(0) - Theta2(0) = 0,
(D(u1))(0) - Z1*(D(u2))(0) = 0,
(D(Theta1))(0) - Z4*(D(Theta2))(0) = 0

For arbitrary values of au, aT, bu, bT there is no reason for E to be equal to 0.
So I use Optimization:-NLPSolve to find a set of values for au, aT, bu, bT such that E is equal, or tlose, to 0.

The other method (resolution of two separate problems, one over [-1, 0] and the other over [0, 1], with formal bv at 0 (let's say a, b, c, ..., for instance u1(0)=a,  u2(0)=a, Theta1(0)=b, Theta2(0)=b, ... in order to satisfy the first relation above) is essentially the same.
Indeed there is no reason that for arbitrary values of a, b, c, ..., the odes are verified over [-1, +1].
For instance you could chose a=0, b=0, c=0, ... and get some pattern of an aggregated solution which trivialy verifies the 4 continuity solutions, verifies the odes in each subdomain, but doesn't verify those ode over the whole (-1, +1) domain.
So here again you will have to sweep the domain a, b, c, ... are assumed to belong to, in order to find some correct values for a, b, c, ....
I found this strategy more complex to manage than the one I developped in my worksheet.

To sum pu, ok to "solved first for -1 to 0 and solved for 0 to 1 and merged both graphs in single plot", but what bv do you take at y=0?

@Andiguys 

display(
  plot3d(
    TRC(tau1, delta), tau1 = 0 .. 1, delta = 0 .. 1
    , axis[3] = [tickmarks = 8] # add more tickmarks for the z axis
  ), 
  # add plane z=0
  plot3d(0, tau1 = 0 .. 1, delta = 0 .. 1, style = wireframe, color = "LightGray"), 
  pointplot3d(P, color = red, symbol = solidcircle, symbolsize = 17), 
  textplot3d([P[], "P"], font = [times, 16], align = {ABOVE, LEFT}), 
  orientation = [165, 75, 0], 
  labels = ['tau[1]', 'delta', 'Profit'], 
  view = [default, default, TRC(1, 1) .. -(1/2)*TRC(1, 1)]  # maximum value onthe z axis
)

I believe you really should give a look to the site plot/options which provides answers to your questions.

You asked "Also, how can I display the z-value at point P?" I recommend you to display these informations in the title and not through textplot3d: indeed text are not managed in a satisfactory way (IMO) when one rotates a 3D  plot.

@salim-barzani 

I dont' see what is complex in what I did; this is elementary algebra.
Never mind, I have update my previous worksheet by adding the conditions in order that a__12 = 0 might not be the only solution to the 29 equations.

help-parameter_sand15_2.mw

@salim-barzani 

Provide an updated worksheet instead of a single code snippet I don't know where to insert in your initial code.

@janhardo @Carl Love

Carl is right.

The problem you solve is "Finding the equation of the sphere which passes through the point A=[5, 1, 4] and whose intersections with planes xy, yz, zx intersect at the origin" (and if they have a common intersection this latter is necessary at the origin for the reason Carl mentioned).

Maybe this is what @yangtheary wanted to say but it is not what it said.

@Andiguys   @acer 

Here is an attempt to synthetize the expression @acer provided

Q_s_enhanced_readibility.mw

(for what it's worth)

@Muhammad Usman 

see assign help page.

Example:  general_procedure_2.mw  (read carefully the 2nd example)