I am trying to write code for analytical solution of 1D heat conduction equation in semi-infinite rod. The analytical solution is given by Carslaw and Jaeger 1959 (p305) as
$$ h(x,t) = \Delta H .erfc( \frac{x}{2 \sqrt[]{vt} } ) $$
where x is distance, v is diffusivity (material property) and t is time. I have written following code in MATLAB for to find heat distribution at 40 m from the input with time. The input is a sine wave at x=0.
clear; close all; clc;
x = 40;       % distance from input to observation point
v = 100;      % diffusivity
p = 300;        % number of steps
Ti = 0;       % Initial head
t =  1:p;    % time steps at which calculation is done, lenth(t) should be
equal to length(To)
delH = sin(linspace(0.1,60,p));

h = zeros(1,p);

for m = 1:p
    func = x/(2*sqrt(v*t(m)));
    h(m) = delH(m) * erfc(func);

hold on;
I am not understanding why the output sinewave keeps its amplitude increasing with time? Should it just mimic the input sine wave with reduced amplitude. The yellow line in the image shows result from finite difference and I expect result from Carslaw to mimic finite difference result. As the book of Carslaw is not available online here is the link to the equation I am using from book of Jan Taler on (chap 14, euqtion 15 on page 357). Modeling analytical solution of 1D heat equation in semi-infinite rod.-carslaw_vs_fd.png