I must be missing something, since the OP gushes over the combination of using PlotVector@MapToBasis@VectorField from the VectorCalculus package, for obtaining cartesian coordinate display of the arrows while using spherical coordinate representation for the vector-field.

But that's just what I'd already shown that VectorField from the Student:-VectorCalculus package provides. (I'd already used the offset for the domain of the azimuth, naturally.)

There must be something I'm missing. Maybe a generality thing.

@Melvin Brown Please put your followups as Replies here, instead of in a separate Question.

@Carl Love Ae you sure that vv gave a way to obtain a colorscheme that was not "zgradient"?

By that I mean: the OP expressed what he wanted as what might be expected from the (currently invalid) syntax,

c := ph+th/16:
plot3d(1, th=0..2*Pi, ph=0..Pi, coords=spherical,
        color=[c,colorscheme=[black, red, yellow, white]]);

I interpreted that to describe a scalar-valued color-function c which would be mapped to the color-gradient formed from those four color anchors. Another (also currently invalid) syntax for that might be, say,

c := ph+th/16:
plot3d(1, th=0..2*Pi, ph=0..Pi, coords=spherical,
       colorscheme=["xyzcoloring",
                    unapply(c,[th,ph,z]), [black, red, yellow, white]]);

Here is another way to obtain the first "replacement" plot in my Answer. This uses indets instead of op to pick out the new color data. And this uses overrideoptions instead of subsop to force the replacement.

restart;
with(plots):
c := ph+th/16:
P := plot3d(1, th=0..2*Pi, ph=0..Pi, coords=spherical, color=c):
Q := plot3d(1, th=0..2*Pi, ph=0..Pi, coords=spherical):
display(Q, overrideoptions=true,
        color=indets(surfdata(op(2,indets(P,specfunc(COLOR))[1]),
                              colorscheme= ["zgradient",[black, red, yellow, white]]),
                     specfunc(anything,COLOR))[1]);

And here is yet another way, using densityplot instead of a second plot3d call and surfdata.

restart;
with(plots):
c := ph+th/16:
Q := plot3d(1, th=0..2*Pi, ph=0..Pi, coords=spherical):
display(Q, overrideoptions=true,
        color=indets(densityplot(c,th=0..2*Pi, ph=0..Pi,
                                colorscheme=["zgradient",
                                             [black, red, yellow, white]]),
                     specfunc(anything,COLOR))[1]);

I expect that there may be yet another way, where the image option is used directly in a single plot3d call. I haven't yet had a second coffee, so I may need to switch variables and use DataTranspose instead of FlipDimension. But the idea is,

restart;
with(plots):
c := ph+th/16:
R:=op(2,indets(densityplot(c, ph=0..Pi, th=0..2*Pi,
                           colorscheme=["zgradient",
                                        [black, red, yellow, white]]),
                         specfunc(anything,COLOR))[1]):
plot3d(1, th=0..2*Pi, ph=0..Pi, coords=spherical,
       image=ArrayTools:-FlipDimension(R,1), style=patch);

Lots of fun.

custom_colorscheme.mw

@anton_dys It would also be possible to construct piecewise coloring functions, for just the original plot3d construction of the sphere, using complicated piecewise formula(s). But to me it seems easier to leverage the colorscheme functionality, and replace the COLOR substructure.

There are other ways to obtain a new COLOR substructure for use as a replacement. Eg, making sure to use the same grid sizes,

restart;
c := ph+th/16:
P2 := plot3d(1, th=0..2*Pi, ph=0..Pi, coords=spherical, grid=[49,49]):
P2 := subsindets(P2,specfunc(anything,{COLOR,COLOUR,SHADING}),()->NULL):
U := plot3d(c, th=0..2*Pi, ph=0..Pi, grid=[49,49],
            colorscheme= ["zgradient",[black, red, yellow, white]]):
subsop([1]=op([1,0],P2)(op([1,..],P2),indets(U,specfunc(anything,COLOR))[1]),P2);

You could even use a different coordinate set for producing the colorscheme colors, or different scheme like "xgradient". (I'm just making this example up...)

restart;
c := ph+th/16:
P2 := plot3d(1, th=0..2*Pi, ph=0..Pi, coords=spherical, grid=[49,49]):
# Next step not strictly necessary here, as P2 has not COLOR substructure.
# But more generally P2 could have arisen elsewhere. So remove any COLOR
# or SHADING from it.
P2 := subsindets(P2,specfunc(anything,{COLOR,COLOUR,SHADING}),()->NULL):
U := plot3d(c, th=0..2*Pi, ph=0..Pi, grid=[49,49], coords=spherical,
            colorscheme= ["xgradient",[black, red, yellow, white]]):
subsop([1]=op([1,0],P2)(op([1,..],P2),indets(U,specfunc(anything,COLOR))[1]),P2);

@minhhieuh2003 You only need to state the followup question and insert a link to any new uploaded worksheet. You don't need to inline the contents of every newly uploaded worksheet.

I advise that you stop using global variables. They are counterproductive and make it awkward. Simply have dsds1d_s return a list of the Rows. And have printdata accept a list of Rows. You can then call printdata with a merged list of the different results from separate calls to dsds1d_s .

goi_maple_new_Q_ac.mw

@torabi For multivariate systems the method fsolve uses is Newton's method.

I doubt that re-implementing Newton's method will help.

But, now that the equations are correct you should be able to try other things.

You could try looking for more suitable ranges for the varables.

You could construct procedures, so as to call fsolve with higher working precision while keeping acceptance/convergence tolerance looser.

Or you could try DirectSearch:-SolveEquations.

@adel-00 Why do you keep ignoring the suggestions which work?

@Preben Alsholm What I was trying to convey, in my Answer above, is that I get the following performance using the Zf(t,f) I had.

plot3d(Zf(t,f), f=1..10, t=0..3, grid=[29,49]);     # 9.0 sec

plot3d(Zf(t,f), t=0..3, f=1..10, grid=[49,29]);     # 3.6 minutes

Using the implementation for Q(f,t) which you've given, I see the following, similar, performance distinction (albeit a smaller absolute difference, but a somewhat comparable relative difference).

CodeTools:-Usage( plot3d(Q(f,t), f=1..10, t=0..3, grid=[29,49]) ):
memory used=11.90MiB, alloc change=0 bytes, cpu time=330.00ms, real time=331.00ms, gc time=30.19ms
CodeTools:-Usage( plot3d(Q(f,t), t=0..3, f=1..10, grid=[49,29]) ):
memory used=49.77MiB, alloc change=0 bytes, cpu time=4.82s, real time=4.82s, gc time=29.40ms

The key aspects are that it seems important to:
1) Have the dsolve parameter be the first/inner plotting variable, so that it changes value less often.
2) Don't set that parameter (using a procedure returned by dsolve) if its value is not different/changed.

Yes, I realize that the axis roles can be interchanged by using the parametric plotd3 calling sequence, while retaining the dsolve-parameter as the inner/first plotting variable for efficiency. But I was trying to keep it simple since the OP has trouble following instructions.
 

@adel-00 No, do it as,

   plot3d(Zf(t,f), f=-1..1, t=0..3, grid=[29,49]);

The procedure Zf sets the value for parameter f, and then calls zfun at the value of t.

The CodeTools:-Usage is just there to show the timing.

Download Summation.mw

@minhhieuh2003 Column span relates to how many columns a cell will span over. It doesn't specify a width.

You might be able to force some kind of increase to width, by padding with a long empty string. But it would be poor, in the sense that it would hardly offer the functionality to get a target pixel size.

I'll repeat an earlier suggestion: Consider using the documented commands in DocumentTools:-Layout instead of all these undocumented commands in Typesetting. It also gives the control you've requested for horizontal alignment as well as width. For example,

goi_maple_ac_EC.mw

@FrancescoLR What I'm expecting is a worksheet or text file, that contains all code, including procedures and data and the function calls that will reproduce your error messages when re-executed.

Use the green up-arrow in the Mapleprimes editor to upload.

@torabi As I wrote before, "Look at the entries in PDE5 and PDE[5], instead of suppressing everything with full colons. The equations you're passing to fsolve contain scalar expressions added to rtables."

You wrote the code which constructs PDE[i], and you haven't explained it. So, you ought to fix it.

@torabi Examine the arguments you are passing to fsolve. They contain some non-scalar nonsense.

@Carl Love He is also attempting to do Fortran code-generation, if I am understanding rightly. But we've been told very little about the data. It might be possible to get some improvement for at least the CG part, perhaps by using Grid on a few chunks.