在上一问的基础上，在上一问的流场中放入一个质量为m的刚性球体，此球体半径为r,试求出以球心为坐标原点的流场的速度函数，给出matlab脚本代码

根据物理学的基本原理，我们可以用势函数来描述该流场，假设流场的势函数为 $\phi(x,y,z)$，则速度场的每个分量可以表示为：$v_x = -\frac{\partial \phi}{\partial x}, v_y = -\frac{\partial \phi}{\partial y}, v_z = -\frac{\partial \phi}{\partial z}$。因为球体是刚性的，所以速度场和势函数应该满足如下条件：

1. 在球体外，速度场满足无旋转条件 $\nabla \times \boldsymbol{v} = \boldsymbol{0}$，即速度场可以表示为梯度的形式 $\boldsymbol{v} = -\nabla \phi$。

2. 在球体内，速度场和势函数应该满足无穿透条件，即球体表面的速度等于球体表面的法向导数和速度场的点积为0。

根据上述条件，我们可以写出势函数的表达式：

$$
\phi(x,y,z) = \begin{cases}
\frac{m}{4\pi r^3}\left(\frac{3r^2}{2}-\frac{r^3}{2\sqrt{x^2+y^2+z^2}}\right), & \text{if } r < \sqrt{x^2+y^2+z^2} \\
\frac{m}{2r}\left(1 - \frac{3}{2}\frac{r^2}{(x^2+y^2+z^2)}\right), & \text{if } r \geq \sqrt{x^2+y^2+z^2}
\end{cases}
$$

其中第一行表示球体外的势函数，第二行表示球体内的势函数。

下面是 Matlab 代码实现：

% 球的半径和质量
r = 1;
m = 1;

% 定义势函数
phi = @(x,y,z) ...
    (r < sqrt(x.^2+y.^2+z.^2)).* ...
    m./(4*pi*r^3) .* (3*r^2/2 - r^3./(2*sqrt(x.^2+y.^2+z.^2))) + ...
    (r >= sqrt(x.^2+y.^2+z.^2)).* ...
    m./(2*r) .* (1 - 3/2*r^2./(x.^2+y.^2+z.^2));

% 绘制流场
[x,y,z] = meshgrid(-2:0.2:2);
u = diff(phi(x,y,z),1,2);
v = diff(phi(x,y,z),1,1);
w = diff(phi(x,y,z),1,3);
quiver3(x(:,:,1),y(:,:,1),z(:,:,1),u(:,:,1),v(:,:,1),w(:,:,1));

运行后可以得到球体流场的速度矢量场图像。