假设你是matlab程序员，已知作用激光功率为P=600w,半径为w=1cm的基模高斯激光，已知岩石样品的密度为ρ=2g/cm3,比热容为C=0.75J/(g.K),热传导系数为K=4.4W/(m.K),假设岩石对光吸收率为η=0.6，岩石为长10cm，宽10cm，高15cm的长方体体,初始温度T0＝300K，边界条件：高度方向两表面为第二类边界条件，根据matlab软件利用有限差分法计算激光以v=0.13cm/s的速度，照射3s后的岩石表面沿移动方向温度场和应力场

题目中给出的参数如下：

- 激光功率 P = 600 W
- 激光半径 w = 1 cm
- 岩石密度 ρ = 2 g/cm3
- 岩石比热容 C = 0.75 J/(g.K)
- 岩石热传导系数 K = 4.4 W/(m.K)
- 岩石光吸收率 η = 0.6
- 岩石长 L = 10 cm
- 岩石宽 W = 10 cm
- 岩石高 H = 15 cm
- 初始温度 T0 = 300 K
- 移动速度 v = 0.13 cm/s
- 照射时间 t = 3 s

我们可以采用有限差分法来计算激光照射岩石后的温度场和应力场。有限差分法是一种数值计算的方法，它将连续的微分方程转化为离散的差分方程，从而得到数值解。

首先，我们需要定义计算区域的大小和离散化步长。由于岩石的尺寸为 10 cm × 10 cm × 15 cm，我们可以将计算区域设置为一个 100 × 100 × 150 的立方体，然后将其离散化为 100 × 100 × 150 个小立方体。为了简化计算，我们可以将激光垂直于岩石的方向设置为 z 轴正方向。

接下来，我们可以使用以下代码来初始化计算区域和参数：

% 初始化计算区域和参数
L = 10;        % 岩石长
W = 10;        % 岩石宽
H = 15;        % 岩石高
dx = 0.1;      % 离散化步长
dt = 0.01;     % 时间步长
nx = L/dx + 1; % x 方向上的离散化点数
ny = W/dx + 1; % y 方向上的离散化点数
nz = H/dx + 1; % z 方向上的离散化点数
T = ones(nx,ny,nz)*300; % 温度场
K = 4.4;       % 热传导系数
C = 0.75;      % 比热容
rho = 2;       % 密度
P = 600;       % 激光功率
w = 0.01;      % 激光半径
eta = 0.6;     % 光吸收率
v = 0.13;      % 移动速度
t = 3;         % 照射时间

接下来，我们可以编写一个函数来计算每个时间步长内的温度场变化。我们可以使用迭代法来解决方程，其中包含激光对岩石的热源项、热传导项和边界条件。

function T = calc_temp(T, dx, dt, K, C, rho, P, w, eta, v, t)
% 计算每个时间步长内的温度场变化
nx = size(T,1);
ny = size(T,2);
nz = size(T,3);
for i = 2:nx-1
    for j = 2:ny-1
        for k = 2:nz-1
            % 计算热传导项
            dTdx2 = (T(i+1,j,k) - 2*T(i,j,k) + T(i-1,j,k))/(dx^2);
            dTdy2 = (T(i,j+1,k) - 2*T(i,j,k) + T(i,j-1,k))/(dx^2);
            dTdz2 = (T(i,j,k+1) - 2*T(i,j,k) + T(i,j,k-1))/(dx^2);
            Q = K * (dTdx2 + dTdy2 + dTdz2);
            % 计算热源项
            G = eta * P / (pi * w^2) * exp(-2*(x^2+y^2)/(w^2)) / rho / C;
            % 计算温度场变化
            T(i,j,k) = T(i,j,k) + dt * (Q + G);
        end
    end
end
% 处理边界条件
T(1,:,:) = T(2,:,:);
T(nx,:,:) = T(nx-1,:,:);
T(:,1,:) = T(:,2,:);
T(:,ny,:) = T(:,ny-1,:);
T(:,:,1) = T(:,:,2);
T(:,:,nz) = T(:,:,nz-1);
end

最后，我们可以使用以下代码来计算激光照射岩石后的温度场和应力场：

% 计算温度场
for i = 1:round(t/dt)
    T = calc_temp(T, dx, dt, K, C, rho, P, w, eta, v, t);
end

% 计算应力场
Sx = zeros(nx,ny,nz);
Sy = zeros(nx,ny,nz);
Sz = zeros(nx,ny,nz);
for i = 2:nx-1
    for j = 2:ny-1
        for k = 2:nz-1
            % 计算应力场
            dTdx = (T(i+1,j,k) - T(i-1,j,k))/(2*dx);
            dTdy = (T(i,j+1,k) - T(i,j-1,k))/(2*dx);
            dTdz = (T(i,j,k+1) - T(i,j,k-1))/(2*dx);
            Sx(i,j,k) = -K * dTdx;
            Sy(i,j,k) = -K * dTdy;
            Sz(i,j,k) = -K * dTdz;
        end
    end
end

温度场和应力场的结果将存储在 T、Sx、Sy 和 Sz 中。您可以使用 MATLAB 中的相关函数来可视化结果。