Hi @문기
I didn't modify your equations (except for removing the 'simplify' part), but you need to provide values for the parameters in the 'ode' function. Some division terms encounter a division-by-zero singularity (referred to as Infinity ∞ in layman's terms) when 
and 
. Therefore, I adjusted the code by using small positive values for tspan(1) and Z0.
and . Therefore, I adjusted the code by using small positive values for tspan(1) and Z0.%% Define the differential equation in a function object
function dZdt = ode(t, Z)
% parameters
theta_m = 1;
alpha = 1;
rho = 1;
del_H = 1;
k = 1;
T1 = 1;
T0 = 0;
% terms
a = theta_m / erf(Z) / (2*sqrt(alpha*t));
c = (1 - theta_m)/(erfc(Z)/(2*sqrt(alpha * t)));
term1 = a - c;
term2 = (sqrt(pi)*rho*del_H/(2*k*(T1-T0))) * exp((Z)^2/(4*alpha^2*t^2));
% ODE
dZdt = term1/term2;
end
tspan = [0.0001 5];
Z0 = 0.0001;
[t, Z] = ode45(@ode, tspan, Z0);
plot(t, Z), grid on, xlabel('t'), ylabel('Z(t)'), title('Time response of Z')
