Is this integral equation valid if applied across a boundary of a magnetized material,(that is one end of the integral is inside the magnetised body and the other is out of it)? Or does it inherit any problems from its differential form?
There is an assumption regarding this boundary condition that is often overlooked for various reasons, but can sometimes become problematic. The limitations of this assumption are discussed in detail in Sections I.5 and I.6 of Classical Electrodynamics, Third Edition by John D. Jackson (i.e., blue cover version). So we start with:
$$
\oint_{C} \ \mathbf{H} \cdot d\mathbf{l} = \int_{S'} \ da \ \left[ \mathbf{j} + \partial_{t} \mathbf{D} \right] \cdot \mathbf{n}' \tag{0}
$$
where $\mathbf{H}$ is the magnetic field (technically, $\mathbf{B}$ is the magnetic induction), $\mathbf{D}$ is the electric displacement, $\mathbf{j}$ is the current density (specifically the macroscopic average current density, see pages 248--258 in Jackson's book for definition and derivation), $S'$ is a closed surface with an outward unit normal $\mathbf{n}'$, and $\partial_{t} = \tfrac{ \partial }{ \partial t }$.
Typically one sweeps the following under the rug, as it were. There can be a surface current density, $\mathbf{K}$, that exists in a thin layer no thicker than one electron skin depth on the surface of the conducting material either caused by time-varying fields or merely present due to some source. In such scenarios, the right-hand side of Equation 0 changes to:
$$
\int_{S'} \ da \ \left[ \mathbf{j} + \partial_{t} \mathbf{D} \right] \cdot \mathbf{t} = \mathbf{K} \cdot \mathbf{t} \ \Delta l \tag{1}
$$
where $\mathbf{t}$ is the unit vector transverse to the surface $S'$ and $\Delta l$ is the scale length of the pill box transverse to the surface $S'$.
There is another issue that arises at such a boundary involving what Jackson refers to as truly microscopic surface charge densities, i.e., $\rho\left( x \right) = \sigma \delta\left( x \right)$ where $\sigma$ is the average surface charge density and $\delta\left( x \right)$ is the Dirac delta function. The idealized scenario we are taught in class is that $\sigma$ exists on the surface and that it has zero thickness, i.e., no charge density inside conductors, only on the surface. However, the truth is that $\rho\left( x \right)$ is confined to within $\pm 2$ angstroms of the ''surface'' of the ionic distribution. For nearly all purposes, there is a discontinuity in the electric field at this boundary but in reality it likely varies over a finite length comparable to a few atomic widths or so.
However the curl is not defined on the boundary so is the loop integral well defined? Another way of saying is that: can I use stokes theorem when the curl of function isn't defined on some region?
It depends on what you are doing. Are you testing what Jackson would call the microscopic or the macroscopic? If the latter, things are much easier and the right-hand side of Equation 0 is can vanish if there is no significant surface current density, $\mathbf{K}$, at the boundary. Even if present, one can still work with Equation 0 or 1, depending on scenario, and get meaningful results for macroscopic approximations.
Side Note
I did not rederive Jackson's definition of microscopic vs macroscopic because it's 10 pages in his book and not really necessary here. It basically involves noting the difference between spatial and temporal ansemble averages and why spatial are the correct choice then lots of details as to why XYZ is okay under WUV limits. The distinction in regards to this question is whether the OP wants to properly model $\mathbf{H}$ across the boundary on scales down to the atomic or if they are okay with the typical, larger scale approximations of micrometers and up (give or take).