Untitled

% Define parameters
rho = 2700; % kg/m^3
Cp = 903; % J/(kg*K)
k = 204; % W/(m*K)
L = 1; % m
T_inf = 298; % K
h = 25; % W/(m^2*K)
q_flux = 1000; % W/m^2
dt = 0.1; % s, time step
dx = 0.01; % m, spatial step size

% Compute number of spatial positions and initialize temperatures
N = round(L / dx);
T = T_inf * ones(1, N); % Initial temperatures

% Compute heat fluxes at boundaries
q_left = q_flux;
q_right = h * (T(end) - T_inf);

% Time loop
for n = 1:100 % Number of time steps
    % Store old temperatures for energy balance check
    T_old = T;

    % Update temperatures
    for i = 2:N-1
        T(i) = T(i) + dt * k / (dx^2) * (T(i-1) - 2*T(i) + T(i+1));
    end

    % Update boundary temperatures
    T(1) = T(1) + dt * q_left / k;
    T(end) = T(end) + dt * q_right / k;

    % Compute energy balance
    energy_old = sum(T_old);
    energy_new = sum(T);
    energy_in = dt * q_left * dx / Cp;
    energy_out = dt * q_right * dx / Cp;
    balance = energy_new - energy_old - energy_in + energy_out;

    % Display temperatures at two time steps
    if n == 50 || n == 100
        disp(['Temperatures at time step ', num2str(n), ':'])
        disp(T)
    end

    % Check energy balance
    if abs(balance) > 1e-6
        disp(['Energy balance check failed at time step ', num2str(n)])
        break
    end
end