@Carl Love Thank you for your help

@Markiyan Hirnyk and Carl, thank you for your help. I got it.

One more question, is there exist a universal method to judge the numerical solution of ode is right?

@rlopez I guess the elliptic problem can be solved like the shoot method with bvp in ode.

That means convert bvp to ivp, and use optimization method to search the solution meet the boundary value.

Of course, the pde is more complicate. It will assume the unknown function on the boundary.

And then minimize the residuals too.

If the solution on the boundary is smooth, then this method will take effect.

@phil76600 I have asked two similar questions.

Unfortunately, the final one model is still unresolved.

http://www.mapleprimes.com/questions/202662-Dsolve-With-Some-Difficult-Ode

http://www.mapleprimes.com/questions/145201-Ode-With-Events

If you got the last model resolved, post it in this forum.

good luck!

@Alejandro Jakubi Thank you very much, my friend.

I'll try it.

If I have problem, I'll ask a new question about that.

thank you for your help.

@Alejandro Jakubi 

One more question.

I found the filedplot3d command could part fill my demand.

for example, plots:-fieldplot3d([f(x,y,z),0,0], x=-1..1,y=-1..1,z=-1..1, color=??)

but how to make the color function?

for example, the min(f(x,y,z)) and max(f(x,y,z)) are red and black. other middle value in the transiton from red to black.

@Alejandro Jakubi thank you very much for  your help.

I'll try it.

@Markiyan Hirnyk thank you for your help.

I'll try it.

@Carl Love thank you for your advice.

In fact, the scalar field which I want to plot is a time dependent field. that means f(x,y,z,t).

I haven't seen an animation about CAT or MRI. But I can image that: an animation of 2d plot.

that's slice of 3d image and can't give a whole figure.

thank you for your suggestion.

@Alejandro Jakubi 

thank you for your reply.

The 3-D (cubic) array of colored points is a good idea.

would you like to give an example about how to color the points array?

thank you in advance. I'll have a try.

@Markiyan Hirnyk thank you for your reply.

But I think you don't understand my question.

A function f(x,y,z) is not a surface in 3d. Am I right?

for example, a density field rho(x,y,z) in a bulk steel is obviously not a surface in 3d.

The implicitplot3d command can't meet my demand.

thank you all the same.

@Carl Love 

Thank you for point that. I'll take care.

@Carl Love 

Thank you very much for your help

@Carl Love I tried remove with match or typematch, but failed.

I use your method to these terms. you're right, I just want to make these derivatives with respect to x equal 0.

here is some approaches from you.

restart;
diff(f1(x,y),x)+diff(f1(x,y),y)+diff(f2(x,y),x)+diff(f2(x,y),y);
remove(has,%,indets(%,'diff'(anything,identical(x))));

@Carl Love thank you.

for problem 1, I want to remove these derivatives with respect to x, I tried remove with match or typematch, but failed.

problem 2, that's exact I want.

thank you for your help