假设你是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后的岩石表面沿移动方向温度场和应力场

首先，需要根据激光功率和激光半径计算出激光强度I，并计算出岩石表面的吸收功率Pabs。根据题目给出的条件，可以得到：
I = P / (π * w^2) = 1.91e7 W/m^2
Pabs = η * I = 1.15e7 W/m^2

接下来，需要利用有限差分法来求解温度和应力场。假设岩石的密度和比热容是均匀分布的，则可以根据热传导方程和连续性方程得到以下计算公式：

dT/dt = K * (d^2T/dx^2 + d^2T/dy^2 + d^2T/dz^2) / (ρ * C)
dSxx/dx = 2 * K * (d^2T/dx^2 - 1/3 * (d^2T/dy^2 + d^2T/dz^2))
dSyy/dy = 2 * K * (d^2T/dy^2 - 1/3 * (d^2T/dx^2 + d^2T/dz^2))
dSzz/dz = 2 * K * (d^2T/dz^2 - 1/3 * (d^2T/dx^2 + d^2T/dy^2))

其中，dT/dt表示温度随时间的变化率，dSxx/dx、dSyy/dy、dSzz/dz分别表示应力场在x、y、z方向上的变化率，d^2T/dx^2、d^2T/dy^2、d^2T/dz^2表示温度在x、y、z方向上的二阶导数。

为了使用有限差分法进行求解，需要将岩石体划分为网格，根据题目给出的条件，可以将岩石体分为100个网格，每个网格的尺寸为0.1cm x 0.1cm x 0.1cm。可以通过以下代码生成网格：

L = 10; % 岩石体长
W = 10; % 岩石体宽
H = 15; % 岩石体高
dx = 0.1; % 网格尺寸
dy = 0.1;
dz = 0.1;
nx = L / dx; % 网格数
ny = W / dy;
nz = H / dz;
x = linspace(0, L, nx);
y = linspace(0, W, ny);
z = linspace(0, H, nz);
[X, Y, Z] = meshgrid(x, y, z);

接下来，需要设置初始温度场和边界条件。根据题目给出的条件，可以将初始温度设为300K，边界条件为第二类边界条件，即表面的温度梯度为0。可以通过以下代码设置初始温度场和边界条件：

T = ones(nx, ny, nz) * 300; % 初始温度场
T(:, :, 1) = 300 + v * Pabs / (K * rho * C); % 底部边界条件
T(:, :, end) = 300; % 顶部边界条件
T(1, :, :) = T(2, :, :); % 左侧边界条件
T(end, :, :) = T(end-1, :, :); % 右侧边界条件
T(:, 1, :) = T(:, 2, :); % 前侧边界条件
T(:, end, :) = T(:, end-1, :); % 后侧边界条件

接下来，需要使用有限差分法来求解温度和应力场，可以通过以下代码实现：

dt = 0.01; % 时间步长
Nt = 300; % 时间步数
for n = 1:Nt
% 计算温度场
d2Tx = (T(3:end, 2:end-1, 2:end-1) - 2*T(2:end-1, 2:end-1, 2:end-1) + T(1:end-2, 2:end-1, 2:end-1)) / dx^2;
d2Ty = (T(2:end-1, 3:end, 2:end-1) - 2*T(2:end-1, 2:end-1, 2:end-1) + T(2:end-1, 1:end-2, 2:end-1)) / dy^2;
d2Tz = (T(2:end-1, 2:end-1, 3:end) - 2*T(2:end-1, 2:end-1, 2:end-1) + T(2:end-1, 2:end-1, 1:end-2)) / dz^2;
dTdt = K * (d2Tx + d2Ty + d2Tz) / (rho * C);
T(2:end-1, 2:end-1, 2:end-1) = T(2:end-1, 2:end-1, 2:end-1) + dt * dTdt;

% 计算应力场
d2Txx = (T(3:end, 2:end-1, 2:end-1) - 2*T(2:end-1, 2:end-1, 2:end-1) + T(1:end-2, 2:end-1, 2:end-1)) / dx^2;
d2Tyy = (T(2:end-1, 3:end, 2:end-1) - 2*T(2:end-1, 2:end-1, 2:end-1) + T(2:end-1, 1:end-2, 2:end-1)) / dy^2;
d2Tzz = (T(2:end-1, 2:end-1, 3:end) - 2*T(2:end-1, 2:end-1, 2:end-1) + T(2:end-1, 2:end-1, 1:end-2)) / dz^2;
Sxx = 2 * K * (d2Txx - 1/3 * (d2Tyy + d2Tzz));
Syy = 2 * K * (d2Tyy - 1/3 * (d2Txx + d2Tzz));
Szz = 2 * K * (d2Tzz - 1/3 * (d2Txx + d2Tyy));
end

最后，可以使用matlab的可视化工具来展示温度和应力场的分布情况。以下是一个简单的示例代码：

figure;
imagesc(squeeze(T(:, :, end/2)));
colorbar;
title('Temperature Distribution');

figure;
quiver3(X, Y, Z, squeeze(Sxx), squeeze(Syy), squeeze(Szz));
title('Stress Field');