I have written a code to reproduce fig 2 in the attached paper but I don't seem to get it right.

I am trying a code to reproduce fig 2 in the attached paper, but I don't seem to get it right. Any help would be apreciated. Thanks

% Material properties
rho = 1.35e3; % Density in kg/m^3
Cp = @(T) (-0.07827 + 0.00223*T - 8.66702e-7*T.^2 + 9.63112e-11*T.^3) * 1e3; % Specific heat in J/g°C converted to J/kgK
alpha = 3.5e3; % Absorption coefficient in m^-1
R = 0.1; % Reflectivity
pulseWidth = 6e-9; % Pulse width in seconds
fluence = 790e-3; % Laser fluence in J/cm^2
I0 = fluence / pulseWidth; % Peak laser intensity

% Spatial domain
L = 1e-3; % Length of the sample in meters (1 mm)
num_elements = 100; % Number of spatial elements
x = linspace(0, L, num_elements); % Spatial grid

% Time domain
total_time = 50e-9; % Total time in seconds
num_time_steps = 500; % Number of time steps
t = linspace(0, total_time, num_time_steps); % Time vector

% Initial temperature distribution
T0 = 20; % Initial temperature in °C
T = T0 * ones(num_elements, 1); % Initial temperature vector

% Thermal diffusivity (D = k / (rho * Cp))
D = 0.4e-4; % Thermal diffusivity in m^2/s
k = D * rho * Cp(T0); % Thermal conductivity in W/m°C

% FEM setup
dx = L / (num_elements - 1); % Spatial step size
dt = total_time / num_time_steps; % Time step size

% FEM matrices
A = zeros(num_elements, num_elements);
B = zeros(num_elements, 1);

% Assembly of FEM matrices
for i = 2:num_elements-1
    A(i, i-1) = k / dx^2;
    A(i, i) = -2 * k / dx^2;
    A(i, i+1) = k / dx^2;
end

% Boundary conditions
A(1, 1) = 1;
A(1, 2) = 0;
A(end, end) = 1;
A(end, end-1) = 0;

% Time-stepping solution
T_history = zeros(num_elements, num_time_steps);
T_history(:, 1) = T;

for n = 1:num_time_steps-1
    B(2:end-1) = alpha * I0 * exp(-alpha * x(2:end-1)) * (1 - R);
    B(end) = T0; % Bottom boundary condition (fixed temperature)
    
    T = T + dt * (A * T + B);
    T(1) = T0; % surface boundary condition 
    T_history(:, n+1) = T;
end

% Plotting the results
figure;
hold on;
plot(t, T_history(1, :), 'r', 'DisplayName', 'Surface');
plot(t, T_history(find(x >= 1e-6, 1), :), 'g', 'DisplayName', '1 \mum');
plot(t, T_history(find(x >= 5e-6, 1), :), 'b', 'DisplayName', '5 \mum');
xlabel('Time (s)');
ylabel('Temperature (°C)');
legend;
title('Temperature profiles at different depths');
hold off;

[Dyuman Joshi]

Please provide information and context as to what the figure is and what formulae/methods are being used to obtain the figure.

The paper does not have any required details about your question or your code, it seems.

[Emmanuel Sarpong]

The paper uses the 1D differential heat conduction equation that's equation 1 and the parameter, I(X,t) is given by equation 2. They apply the finite element method (FEM) to do the computations. The quantities, density (P) in g/cm^3, heat capacity(Cp) in Jg^-1K^-1, diffusivity (D) in cm^2/s, absorption coefficient (alpha) at 532 nm given in cm^-1. The thermal conductivity is the product of diffusivity, density and the specific heat.