# Underflow and Machine Precision

I am using Mathematica to solve a system of differential equations and evaluate and plot the solutions on an interval. The solutions have non-zero solutions, but Mathematica underflows in a portion of the expression and as a result provides an incorrect zero result. It seems like there is a way to overcome this but I seem unable to implement it in my code. My code is as follows:

Clear["Global*"]

Ac = 1.9793*10^(-6); (*m^2*)
k = 7.1;(*W/m*K*)
P = 0.005;(*m*)
U = \
0.04;(*m/s*)
Cp = 553;(*J/kg*K*)
q = 300/(P*h);(*W/m^2*)
\[Sigma] =
5.6703*10^(-8);(*m^2*kg/(s^2*K^1)*)
\[Epsilon] = \
0.35;(*dimensionless*)
To = 295;(*K*)
Ti = 295;(*K*)
Tl = 2273;(*K*)
\
h = 0.025;(*m*)
l = 0.05;(*m*)
\[Rho] = 4506;(*kg/m^3*)

eqns = {(Ac*k)*x''[t] - (\[Rho]*Ac*U*Cp)*x'[t] + P*q == 0,
(Ac*k)*y''[t] - (\[Rho]*Ac*U*Cp)*y'[t] - hr[Ta]*y[t] == 0,
x[0] == To, y[l] == Tl, x[h] == y[h], x'[h] == y'[h]};
s = DSolve[eqns, {x[t], y[t]}, t][[1]];
Y[t_] = y[t] /. s // FullSimplify
X[t_] = x[t] /. s // FullSimplify

hr[Ta] = 300; NotDone = True; n = 1;
While [NotDone, {Y[t]; Ta[Y] = (1/(l - h))*Integrate[Y[t], {t, h, l}];
hrn[Ta] = \[Sigma]*\[Epsilon]*P*(Ta[Y] + Ti) (Ta[Y]^2 + Ti^2);
err = Abs[(hr[Ta] - hrn[Ta])/hr[Ta]]; NotDone = err > 0.001;
Print[Ta[Y]];
hr[Ta] = hrn[Ta];  n++}]
`

The average value of the function returns a correct result of approx. 1740 but evaluating the function directly returns a zero value and the error "Exp[-877.895] is too small to represent as a normalized machine number; precision may be lost."

Is there a simple way to circumvent this?