怎样用matlab画一个二重积分（含两个未知数）关于f和两个未知数的三维图像

怎样用matlab画一个二重积分（含两个未知数）关于F和两个未知数的三维图像？

r和θ是未知参数 就是画出F、r、θ三者的三维图像？
我参照网上的例子写了一个，但总报错？ 谁能帮我追加100分。
clear all; close all
syms n;
i = 1;
r = 0.01;
g = 0;
es = [0:0.001:0.01].'; % 变量e的采样序列
hs = [0:0.001:0.02].'; % 变量h的采样序列
[ES, HS] = meshgrid(es, hs); % 生成平面上的小矩形顶点坐标值的矩阵
YS=AS=BS=CS=DS=ES=FS=AA=A=BB=CC=DD=EE=FF=0*ES;
for g = 1:length(es) %length(es）es长度
e = es(g);
for k = 1:length(hs) %length(hs）hs长度
h = hs(k);
f= @(x,y)(1/((e*cosh-x*cosy)^2+(e*sinh-x*siny)^2+0.005^2)-1/((e*cosh-x*cosy)^2+(e*sinh-x*siny)^2+0.025^2))*(e*cosh-x*cosy)*x; %x=r1 y=o1 e=r h=0
AS(k,g) = -symsum(quad2d(f,n*pi/3-pi/12,n*pi/3+pi/12,0.09,0.14),0,5); %前两项是外层上下限
BS(k,g) = -symsum(quad2d(f,n*pi/3+pi/12,n*pi/3+pi/4,0.09,0.14),0,5); %前两项是外层上下限
YS(k, g) = AS(k,g)+BS(k,g);
end
end
YS = 8.2*e^4*YS;
surf(ES,HS,YS);
求和符号下面打错了，n从0到5

楼上两位的回答都很精彩，但也都存在一些问题：

1、chinasunsunsun的回答主要有以下问题（查了好半天才把问题查全）：

（1）楼主给的公式有两项积分，但您只计算了一项，而且是把第一项的被积函数与第二项的积分限混在一起使用；

（2）楼主的n取值是0~5，您只取了0~4；

（3）画图的方式：很明显θ作为角度量，画图应该按照柱面坐标，但您是按照直角坐标画的（相应的，θ取值范围也不恰当）。

2、AliKsai回答的内容非常丰富，但楼主想要的是F对r、θ的三维图像而已，您前面的大量讨论都不能反映楼主的真正需求，最后的那张图似乎有那么点意思了，但终究还是难见全貌，正确性也可疑，而且过程过于繁琐。

以上点评纯粹出于技术探讨，希望不要介意。

我写的代码如下，如果有问题欢迎批评：

f1=@(t,r,t1,r1)(-1./((r*cos(t)-r1.*cos(t1)).^2+(r*sin(t)-r1.*sin(t1)).^2+0.005^2).^(3/2)+...
    1./((r*cos(t)-r1.*cos(t1)).^2+(r*sin(t)-r1.*sin(t1)).^2+0.025^2).^(3/2)).*(r*cos(t)-...
    r1.*cos(t1)).*r1;
f2=@(t,r,t1,r1)(1./((r*cos(t)-r1.*cos(t1)).^2+(r*sin(t)-r1.*sin(t1)).^2+0.005^2).^(3/2)-...
    1./((r*cos(t)-r1.*cos(t1)).^2+(r*sin(t)-r1.*sin(t1)).^2+0.025^2).^(3/2)).*(r*cos(t)-...
    r1.*cos(t1)).*r1;
An = @(n,t,r) 8.2e4*dblquad(@(r1,t1)f1(t,r,t1,r1),0.09,0.14,(n/3-1/12)*pi,(n/3+1/12)*pi) ...
    + 8.2e4*dblquad(@(r1,t1)f2(t,r,t1,r1),0.09,0.14,(n/3+1/12)*pi,(n/3+1/4)*pi);
F = @(t,r)sum(arrayfun(@(n)An(n,t,r),0:5));

r = linspace(0,.1,20);
t = linspace(0,2*pi,20);
[t,r]=meshgrid(t,r);
z = arrayfun(F,t,r);
[x,y,z]=pol2cart(t,r,z);
surf(x,y,z)

其实我算的是如上图永磁体盘12块两种极性磁铁对导磁体盘产生的磁场强度，我把解析式子推导出来了，平时都用有限元仿真软件做，没解析做过，那两个式子分别是永磁体上端面和下端面对A点产生的磁场强度。r和θ变化相当于导磁体盘底部A点的变化，我需要画出A点在不同位置时（r和θ变化）的磁场强度。您的回答很专业，待我验证后给分，谢谢您的慷慨无私的回答！

谢谢你对问题的背景做了说明。我回答过的问题当中，有不少不知道到底是什么应用背景，只是从数学或编程的角度给出求解方法。能对背景稍微有点了解会感觉具体了很多。

我现在取的点偏少，可以修改

    r = linspace(0,.1,20);
    t = linspace(0,2*pi,20);
为

    r = linspace(0,.1,30);
    t = linspace(0,2*pi,60);
得到的结果更平滑一些（当然，也可以取更多点，但会大大增加计算量）。

另外，对于磁场强度的表现，也可以考虑用伪彩图（pcolor），它是以色彩表示强度：

pcolor(x,y,z)
axis square
colorbar
shading interp

还有，你现在要表现的磁场强度只是一个标量，需要的话，可以同时计算出各处的磁力线方向，进而用三维的场图去表现。

您好，您的回答对我启迪很大，还要麻烦您看下下面的式子。和我按您的程序改的，画出的图有点奇怪，您看下我写的是否正确，还请不吝赐教，谢谢！由于字数限制，我只好截图。

本来想请把你现在写的代码发给我（如果这里受字数限制，可通过网盘或私信），一来避免我重新录入，二来看看你现在的图到底什么样，有问题的话也好查原因。但等了很久也没等到，所以还是自己改了一下代码。

如果按照我原来的取值范围，即r=0~0.1，θ=0~2*pi，结果如下（看上去似乎和你的磁铁布局更吻合了）：

如果按照你代码中的取值范围，即r=0~0.2，θ=0~pi/3，结果如下：

你画出的图是这样的吗？具体的物理概念我没仔细想，但看上去似乎没什么太奇怪的啊？