  We study the equilibrium configurations of relativistic Neutron Stars(NS)
with a polytropic model in a $f(R,T)=R+2\lambda T+\xi T^{2}$ gravity.We
investigate the neutron star properties and their dependence on $\lambda$ and
$\xi$ corresponding to different central densities ($\rho_c$) of the NS. For
$\lambda = 0,-1,-3,-5$ with $\xi=0$ and
$\rho_c=1.5\times10^{18}~\rm{kg~m^{-3}}$, we find the maximum mass of the NS as
$M = 1.06 M_\odot$, $1.19 M_\odot$ $1.61 M_\odot$ and $2.47~M_\odot$
corresponding to the radius($R$) $10.409$ km, $10.737$ km, $11.461$ km and
$12.119$ km. This higher value of NS mass can be compared with gravitational
wave data(GW170817). For given $\lambda =-6$ and $\xi = 0$, we find that as
$\rho_c$ increases from $\rho_c=1.1 to 1.6 \times 10^{18}~\rm{kg~m^{-3}}$, the
maximum mass of the NS decreases from $4.19 M_\odot$ to $3.23 M_\odot$ while
it's radius $R$ decreases $13.86 \rm{km}$ to $11.54 \rm{km}$. With the fixed
value of $\xi = 10^{-27}$ and $\lambda = 0,-1,-3,-5$, we find the maximum mass
$M =1.06 M_\odot$,$1.34 M_\odot$,$1.89 M_\odot$ and $3.39~M_\odot$
corresponding to the radius $R = 10.409$ km, $10.843$ km, $11.549$ km and
$11.680$ km. respectively. Taking our observational constraints i.e. GW170817
(BNS Merger) mass - radius data, observed pulsars PSRJ1614-2230, PSRJ0348+0432
maximum mass - radius data; we found that posterior distribution plot of mass
$\&$ radius gives good result and the corner plot of modified gravity
parameters $\lambda$ and $\xi$ are giving very good posterior results. So, for
a range of values of $\lambda$ with $\xi=0 (\neq 0)$, we found that the mass
$M$ and the radius $R$ of the NS lie within the range given by the GW170817
gravitational wave data given by LIGO, Pulsars $\&$ Millisecond Pulsars data
and the NICER (Neutron star Interior Composition ExploreR) mass-radius data
given by NASA.
