function T = efdiff(T0,dt,tmax,dz,zmax)
%
% EFDIFF
%
%    T = efdiff(T0,dt,tmax,dz,zmax)
%    Explicit formulation of finite difference simulation for
%    the heat conduction problem.  T0 is the uniform initial
%    temperature of the tissue. dt is the time step, tmax is 
%    the final time value, dz is the distance increment, zmax
%    is the maximum depth of simulation.  Other parameters
%    describing the laser source, and the tissue properties
%    must be changed in the beginning of the m-file itself.
%
%                                                     RJM 11/28/94

mu = .5;  %    (1/cm)
%     rho = 1000; %  (1000 kg/m^3)
%     c = 4180; %    (J/°K/s/kg)
%     rhoc=rho*c;
rhoc = 4e+6;    %    (J/m^3)
k = .6; %      (W/m/°K)
I0 = 1.4*(100^2); %    (W/m^2) =  (100^2)*(W/cm^2)
alpha = k/rhoc;


S = dz^2/dt/alpha    % print stability criterion

a = floor(tmax/dt);
b = floor(zmax/dz);
z = dz:dz:zmax;

T=zeros(a,b);
T(1,:) = T0*ones(1,b);  % I.C.

for k=1:a,

 Tleft =  [T(k,2:b)   T0];      % left shift   (constant T at zmax) 
 Tright = [T(k,1) T(k,1:b-1)];  % right shift  (surface not constant)
 Told =   T(k,:);               % unshifted
 Heatin = mu*(100)*I0*exp(-mu*100*z);

 T(k+1,:) = Told + alpha*dt/(dz^2)*(Tleft-2*Told+Tright) + dt/rhoc*Heatin;

end