已知岩石样品的密度为ρ=2g/cm3,比热容为C=0.75,热传导系数为K=4.4,假设岩石对光吸收率为η=0.6,初始温度T0＝300K．利用matlab根据拉普拉斯求沿x轴速度v移动的基模高斯激光辐照岩石温度场，再根据热位移平衡方程求得应力场，请给出代码。

以下是基于Matlab的有限元计算代码，可求解该问题的温度场和应力场。

clc; clear; close all;

%% 岩石样品参数
rho = 2;        % 密度 (g/cm^3)
C = 0.75;       % 比热容 (J/g.K)
K = 4.4;        % 热传导系数 (W/m.K)
alpha = K / (rho * C);  % 热扩散系数 (m^2/s)

%% 光学参数
P = 600;        % 激光功率 (W)
w = 0.01;       % 激光半径 (m)
eta = 0.6;      % 岩石对光吸收率

%% 岩石样品尺寸
Lx = 0.1;       % 长度 (m)
Ly = 0.1;       % 宽度 (m)
Lz = 0.15;      % 高度 (m)
Nx = 101;       % 离散网格数
Ny = 101;
Nz = 151;
dx = Lx / (Nx - 1);
dy = Ly / (Ny - 1);
dz = Lz / (Nz - 1);
x = linspace(0, Lx, Nx);
y = linspace(0, Ly, Ny);
z = linspace(0, Lz, Nz);

%% 时间参数
tmax = 3;       % 最大时间 (s)
dt = 0.001;     % 时间步长 (s)
Nt = round(tmax / dt) + 1;
t = linspace(0, tmax, Nt);

%% 边界条件
T0 = 300;       % 初始温度 (K)
q_laser = P / (pi * w^2) * eta;    % 激光热流密度 (W/m^2)
T = T0 * ones(Nx, Ny, Nz);
T(:, :, 1) = T0;
T(:, :, end) = T0;
T(:, 1, :) = T0;
T(:, end, :) = T0;
T(1, :, :) = T0;
T(end, :, :) = T0;

sigma_xx = zeros(Nx, Ny, Nz);
sigma_yy = zeros(Nx, Ny, Nz);
sigma_zz = zeros(Nx, Ny, Nz);
tau_xy = zeros(Nx, Ny, Nz);
tau_xz = zeros(Nx, Ny, Nz);
tau_yz = zeros(Nx, Ny, Nz);

%% 求解温度场
for k = 2 : Nt
    % 计算激光热流密度
    q = zeros(Nx, Ny, Nz);
    for i = 1 : Nx
        for j = 1 : Ny
            q(i, j, end) = q_laser;
        end
    end
    
    % 计算温度场
    for i = 2 : Nx-1
        for j = 2 : Ny-1
            for l = 2 : Nz-1
                T(i, j, l) = T(i, j, l) + alpha * dt * (T(i+1, j, l) - 2*T(i, j, l) + T(i-1, j, l)) / dx^2 ...
                    + alpha * dt * (T(i, j+1, l) - 2*T(i, j, l) + T(i, j-1, l)) / dy^2 ...
                    + alpha * dt * (T(i, j, l+1) - 2*T(i, j, l) + T(i, j, l-1)) / dz^2 + q(i, j, l) * dx^2 / (rho * C);
            end
        end
    end
    
    % 边界条件
    T(:, :, 1) = T0;
    T(:, :, end) = T0;
    T(:, 1, :) = T0;
    T(:, end, :) = T0;
    T(1, :, :) = T0;
    T(end, :, :) = T0;
    
    % 可视化温度场
    fig1 = figure(1);
    set(gcf, 'Position', [200, 200, 800, 400]);
    subplot(1, 2, 1);
    [X, Y] = meshgrid(y, x);
    surf(X, Y, T(:, :, end)', 'EdgeColor', 'none');
    xlabel('y (m)');
    ylabel('x (m)');
    zlabel('T (K)');
    title(['T(x,y) at t = ', num2str(t(k)), 's']);
    axis([0 Ly 0 Lx 0 1000]);
    colorbar;
    
    % 计算应力场
    for i = 2 : Nx-1
        for j = 2 : Ny-1
            for l = 2 : Nz-1
                sigma_xx(i, j, l) = -K * (T(i+1, j, l) - T(i-1, j, l)) / (2 * dx);
                sigma_yy(i, j, l) = -K * (T(i, j+1, l) - T(i, j-1, l)) / (2 * dy);
                sigma_zz(i, j, l) = -K * (T(i, j, l+1) - T(i, j, l-1)) / (2 * dz);
                tau_xy(i, j, l) = -K * (T(i+1, j+1, l) - T(i+1, j-1,