function dy =Solve_Two_Temperature(t,y)
%t 代表時間變化量
%y 代表所求的電子及晶格溫度
dy = zeros(2,1);    % a column vector
g=4.2e17;%耦合系數
Ce=194;%電子熱容
Ci=3.6e6;%晶格熱容
R=0.431;%反射率
delta=24e-9;%趨膚深度
tao=70e-15;%脈寬
J=4e-9/(pi*20e-12);%能量密度
A=-g/Ce;%系數
 
S=0.94*(1-R)/(tao*delta)*J*exp(-4*log(2)*(t/tao-2).^2);%體激光源項
%dt = zeros(2,1); 
% dt(1)=0.0;
% dt(2)=700e-15;

% S(1)=0.94*(1-R)/(tao*delta)*J*(exp(-4*log(2)*(t/tao-2).^2)+exp(-4*log(2)*((t-dt(1))/tao-2).^2));%雙脈沖體激光源項
% S(2)=0.94*(1-R)/(tao*delta)*J*(exp(-4*log(2)*(t/tao-2).^2)+exp(-4*log(2)*((t-dt(2))/tao-2).^2));
%dy(1) =A*1/y(1)*(y(1)-y(2))+S(1);%方程1 ；y1 電子溫度；y2 晶格溫度
dy(1) =A*1/y(1)*(y(1)-y(2))+S;
dy(2) =g*(y(1)-y(2))/Ci;%方程2