It depends entirely on what $F$ is. The $F = ma$ equation is not a definition! $F$ is defined separately as some function of the positions and velocities that make up the system. Which function $F$ is determines what the dynamics are. And the symmetries of that dynamics hinge on what $F$ is.
In particular, if $F = 0$, then you have just the law of inertia. The symmetry group for it, with three-dimensional space, is $\text{SL}(5)$; and with some minor fudges, you can get it up to $\text{GL}(5)$. The next simplest case are for velocity-independent forces, particularly where the force $F(x) = A + Bx$ (which is in vector-matrix form, if $x$ is two or more dimensional) is an affine function of the coordinates. But even that is involved.
In contrast, what you ended up showing is what transformation laws $F$ ought to satisfy, if the equation, itself, is to be covariant with respect to a given set of symmetries. That restricts the range of choices for $F$. In fact, that's one of the key ways to lay out a classification for possible force laws. It's a model-building technique for dynamics. This holds true not just for mechanics, but for field theory.
I'll give you another example. The symmetries for Maxwell's equations, when written this way:
$$
∇· = σ, \hspace 1em ∇× + \frac{∂}{∂t} = -, \\
∇· = ρ, \hspace 1em ∇× - \frac{∂}{∂t} = +,
$$
(where I expanded it to include non-zero magnetic charges and currents) is huge! It includes all coordinate transformations for $(t,x,y,z)$ - the full transformation group for 4-dimensional manifolds. And, it includes a lot more besides that. Even if you set the magnetic sources to zero, it's still very big. One symmetry - with zero magnetic sources - is $(, ) → ( + θ, - θ)$ (with all other fields remaining unchanged). For lack of a better name, we could call it the "axion transform". It's only when you add in the "constitutive relations"
$$ = ε_0, \hspace 1em = μ_0,$$
that the symmetry group is reduced to something much smaller that includes the Lorentz transforms (as a subgroup), but excludes the Galilean transforms.
The correct way to handle the issue of symmetries is within the Jet Bundle formalism. It is a standard exercise, in it, to determine the symmetries of a given system. That system can involve partial and/or total differentials and/or algebraic relations. For all intents and purposes, even though it is not advertised that way, the Jet Bundle formalism is none other than the Universal Theory Of Differential Equations.
Do you really want me to show you how it pans out for the $F = 0$ case? Do I need to pull Jack Nicholson on you about the truth? Even for that: it's a long derivation! And it needs exposition in a run-up to it, because hardly anybody knows anything about Jet Bundles or has ever seen it actually put to use.
Well, step back and make way, because here it comes...
Outline Of The Jet Bundle Formalism
A system specified by equations or relations involving partial differentials is understood to characterize a trajectory or unfolding over a base space, $M:(x^m)$ expressed locally by coordinates $(x^m)$ that comprise the independent variables of the system. The trajectory or unfolding is described locally by a set of dependent coordinates arranged in a manifold $Q:(q^a)$ called a fibre and referred to as the configuration space. At each $x ∈ M$ is attached a copy of this space $Q_x: (q^a_x)$, the instantaneous configuration space at $x$; and their combination gives us a configuration bundle $Y = \bigcup_{x∈M}\left(\{x\}×Q_x\right)$ with local coordinates $x$. Each configuration or trajectory is then identified locally as a function from the base space to its fibres called a [cross-]section.
$$q: x ∈ M ↦ \left(x, q^a(x)\right) ∈ \{x\}×Q_x ⊆ Y.$$
For convenience, I will introduce the spaces where things live - as above - by listing the coordinates natural to them.
For field theory, there are four independent variables, and $M:(t,x,y,z)$ has 4-dimensions and sections describe both an unfolding within space $(x,y,z)$ and evolution over time $t$ of the fields in question. The fields and their components live in $Q$. The relevant equations are partial differential equations. For mechanics, there is just one independent variable $t$, normally identified as "time" and $M:(t)$ has 1-dimension. Sections describe trajectories and evolution over time and $Q$ contains the full space of objects whose evolution is being laid out. Here: the equations are total differential equations. That's the case of interest here. Mechanics can be regarded as field theory over the "time manifold".
A unified framework may be established for a system of partial differential equations – and for a more general system, such as one which includes inequalities – by “removing the calculus” from the system and converting it into purely algebraic form. This removal is accomplished by treating the derivatives as separate configuration variables. The calculus is recovered by subjecting sections to a restriction that we term here the kinematic constraint, each such constraint embodying what we term its associated kinematic law.
Thus, the first order jet bundle $J^1Y:(x^m, q^a, v^a_m)$ contains extra coordinates $v^a_m$ representing the components of the velocity/gradient $∂_m q^a$. The kinematic constraint is that we associate each section with its jet prolongation:
$$j^1 q: \left(x^m\right) ↦ \left(x^m, q^a, v^a_m\right) = \left(x^m, q^a(x), ∂_mq^a(x)\right).$$
The constraints embody the first order kinematic law: $v^a_m = ∂_mq^a$. The result of separating out this law from the underlying geometry is to turn a system of first order partial differential equations over the bundle $Y$ into a purely algebraic system over the jet bundle $J^1Y$.
A similar extension to $J^1Y: \left(x^m, q^a, v^a_m, a^a_m\right)$ to the second order jet bundle allows us to handle second order prolongations
$$j^2q:\left(x^m\right) ↦ \left(x^m, q^a(x), \frac{∂q^a}{∂x^m}(x), \frac{∂^2q^a}{∂x^m∂x^n}(x)\right),$$
and to treat second order laws of evolution, by adding in extra coordinates $a^a_{mn}$ for the acceleration, which are subject to the kinematic law $∂_mv^a_n = a
^a_{mn}$.
Thus, a second order system
$$F\left(x, q^a(x), \frac{∂q^a}{∂x^m}(x), \frac{∂^2q^a}{∂x^m∂x^n}(x)\right) = 0$$
is reduced to a combination of an algebraic system
$$F(x, q, v, a) = 0$$
in the jet bundle $J^2Y%$, with the “calculus” part of the description
$$\frac{∂q^a}{∂x^m} = v^a_m, \hspace 1em \frac{∂v^a_m}{∂x^n} = a^a_{mn},$$
(which we may think of as the kinematic constraints) pulled back into the infrastructure of the jet bundle formalism itself.
Although we will not go through the details of their development here, a similar set of extensions applies to vector fields over the configuration bundle. Thus, corresponding to a vector field
$$U = U^m \frac{∂}{∂x^m} + U^a \frac{∂}{∂q^a}$$
are its jet prolongations, illustrated here for the first and second order:
$$
\text{pr}^1U = U^m\frac{∂}{∂x^m} + U^a\frac{∂}{∂q^a} + U^a_m\frac{∂}{∂v^a_m} \\
\text{pr}^2U = U^m\frac{∂}{∂x^m} + U^a\frac{∂}{∂q^a} + U^a_m\frac{∂}{∂v^a_m} + U^a_{mn}\frac{∂}{∂a^a_{mn}}, \\
$$
with the extra components $U^a_m$ and $U^a_{mn}$ expressed in terms of $U^m$, $U^a$ and its derivatives. (Footnote: in here and the following, I will use the Einstein summation convention: where a repeated index in a monomial - once in the up position and once in the down position - is understood to be summed over.)
Our task at hand only concerns a system described by ordinary differential equations, where the base space $M:(t)$ reduces to the one-dimensional "time manifold" arena we regard as time and the fibre coordinates $Q:\left(x^μ: μ = 0,⋯,n-1\right)$ yield the extra coordinates for the configuration bundle $Y:\left(t,x^μ\right)$, in which the dynamics of a system described by ordinary differential equations may be formalized. These $x$'s are not spatial coordinates, both the coordinates describing a single object and its trajectory. Correspondingly, for the first and second order jet bundles, we have the coordinates
$$J^1Y: \left(t, x^μ, v^μ\right), \hspace 1em J^2Y: \left(t, x^μ, v^μ, a^μ\right),$$
with $v^μ$ and $a^μ$ respectively representing the components of velocity and acceleration. In mechanics, these are the spaces where velocities and accelerations - the first and second order rates of change of a system's attributes - live. Here: the system is just a single body, and its attributes are the coordinates for its center of mass position.
Jet Bundle Formulation of the Law of Inertia
We may first reformulate the law of inertia in general form for a space of dimension $n ∈ \{1,2,3,⋯\}$ as:
$$\frac{d^2x^μ}{dt^2} = 0.$$
Then, we may pose the question of what its symmetries are. In the framework of second order jet bundles, the equations of motion for the law of inertia reduce to the system
$$a^μ = 0; \hspace 1em \frac{dx^μ}{dt} = v^μ, \hspace 1em \frac{dv^μ}{dt} = a^μ.$$
The symmetry may then be given as the flow of the following vector field
$$δ = δt(x,t) \frac{∂}{∂t} + δx^μ(x,t)\frac{∂}{∂x^μ}.$$
In addition, we can also include $t$ in the set of coordinates by writing $x^0 = t$.
The set of vector fields, that make up the symmetries, form the Lie algebra of the respective symmetry group.
We will seek out the general solution, where no relation is initially assumed between $t$ and the $\left(x^μ\right)$; and then later apply the restriction $x^0 = t$ by setting $δx^0 = δt$. We could go even further, and rewrite the equation as $a^μ = N(t)v^μ$ in order to make the system fully reparametrization invariant. As it is, we will only find a small subset of this invariance group consisting only of rescaling and translations of $t$.
The second order jet prolongation
$$\text{pr}^2δ = δ + δv^μ\frac{∂}{∂v^μ} + δa^μ\frac{∂}{∂a^μ}$$
is given by
$$δv^μ = \frac{D}{Dt}δx^μ - v^μ\frac{D}{Dt}δt, \hspace 1em δa^μ = \frac{D}{Dt}\frac{D}{Dt}δx^μ - 2a^μ \frac{D}{Dt}δt - v^μ\frac{D}{Dt}\frac{D}{Dt}δt,$$
where
$$\frac{D}{Dt} = \frac{∂}{∂t} + v^μ\frac{∂}{∂x^μ} + a^μ\frac{∂}{∂v^μ}.$$
We don’t actually need $δv^μ$ for what follows. As for $δa^μ$, after the reduction
$$\frac{D}{Dt}\frac{D}{Dt} = \frac{∂}{∂t}\frac{∂}{∂t} + 2v^μ∂_μ\frac{∂}{∂t} + v^μv^ν∂_μ∂_ν + a^ν∂_ν$$
we get
$$\begin{align}
δa^μ &= \frac{∂}{∂t}\frac{∂}{∂t}δx^μ + v^ν\left(2∂_νδx^μ + δ^μ_ν\frac{∂}{∂t}δt\right) + v^νv^ρ\left(∂_ν∂_ρδx^μ + \left(δ^μ_ν∂_ρ + δ^μ_ρ∂_ν\right)\frac{∂}{∂t}δt\right) \\
&- v^μv^νv^ρ∂_ν∂_ρδt + a^ν\left(∂_νδx^μ - 2δ^μ_ν\frac{∂}{∂t}δt\right) - a^νv^ρ\left(2δ^μ_ν∂_ρ + δ^μ_ρ∂_ν\right)δt.
\end{align}$$
I'm using the Kronecker delta, defined by $δ^μ_ν = 1$ if $μ = ν$ and $δ^μ_ν = 0$ if $μ ≠ ν$.
The symmetry we seek is then that $δa^μ = 0$, where the equality is taken subject to the constraint $a^μ = 0$ (i.e. "on-shell"). Introducing Lagrange multipliers $Λ^μ_ν$ and $Λ^μ_{νρ}$, we may write:
$$δa^μ = a^ν\left(Λ^μ_ν + Λ^μ_{νρ}v^ρ\right).$$
From this follows the system:
$$
\frac{∂}{∂t}\frac{∂}{∂t}δx^μ = 0, \hspace 1em -∂_ν∂_ρδt = 0, \\
\frac{∂}{∂t}\left(2∂_νδx^μ - δ^μ_ν\frac{∂}{∂t}δt\right) = 0, \hspace 1em ∂_νδx^μ - 2δ^μ_ν\frac{∂}{∂t}δt = Λ^μ_ν, \\
∂_μ∂_νδx^μ - \left(δ^μ_ν∂_ρ + δ^μ_ρ∂_ν\right)\frac{∂}{∂t}δt = 0, \hspace 1em -\left(2δ^μ_ν∂_ρ + δ^μ_ρ∂_ν\right)δt = Λ^μ_{νρ}.
$$
After solving for the Lagrange multipliers
$$Λ^μ_ν = ∂_νδx^μ - 2δ^μ_ν\frac{∂}{∂t}δt, \hspace 1em Λ^μ_{νρ} = -\left(2δ^μ_ν∂_ρ + δ^μ_ρ∂_ν\right)δt,$$
we get
$$
\frac{∂}{∂t}\frac{∂}{∂t}δx^μ = 0, \hspace 1em ∂_ν∂_ρδt = 0, \\
\frac{∂}{∂t}\left(2∂_νδx^μ - δ^μ_ν\frac{∂}{∂t}δt\right) = 0, \hspace 1em ∂_ν∂_ρδx^μ = \left(δ^μ_ν∂_ρ + δ^μ_ρ∂_ν\right)\frac{∂}{∂t}δt.
$$
We conclude from the second of these equations that
$$δt = T(t) + U_μ(t)x^μ.$$
Upon substitution into the last of these equations, we get:
$$∂_ν∂_ρδx^μ = δ^μ_νU_ρ + δ^μ_ρU_ν.$$
From this, we conclude that
$$δx^μ = X^μ(t) + V^μ_ν(t)x^ν + U_ν(t)x^νx^μ.$$
Upon substitution into the remaining equations, this leads to the following:
$$\ddot{X^μ} + \ddot{V^μ_ν}x^ν + \dddot{U_ν}x^μx^ν = 0, \hspace 1em 2\dot{V^μ_ν} + 2x^μ\ddot{U_ν} + δ^μ_ν\ddot{U_ρ}x^ρ - δ^μ_ν\ddot{T} = 0.$$
If we contract the last equation with respect to $μ$ and $ν$, setting $n = δ^μ_μ$, then we get:
$$\frac{n\ddot{T} - 2\dot{V^ρ_ρ}}{n+2} = \ddot{U_ρ}x^ρ.$$
Since $T(t)$, $V^μ_ν(t)$ and $U_μ(t)$ are independent of the $\left(x^μ\right)$, then this means that
$$\ddot{U_ρ} = 0, \hspace 1em \dot{V^ρ_ρ} = \frac{n}{2}\ddot{T}.$$
Upon further substitution, this leads to the following reduced system:
$$\left\{\begin{matrix}
\ddot{X^μ} + \ddot{V^μ_ν}x^ν = 0 \\
2V^μ_ν - δ^μ_ν\ddot{T} = 0
\end{matrix}\right\} ⇒
\left\{\begin{matrix}
\ddot{X^μ} = 0, & \ddot{V^μ_ν} = 0 \\
\ddot{U_ν} = 0, & \dot{V^μ_ν} = \frac{1}{2}δ^μ_ν\ddot{T}
\end{matrix}\right\}.
$$
Thus, we find
$$
X^μ(t) = X^μ_0 + X^μ_1 t, \hspace 1em T(t) = T_0 + T_1 t + T_2 t^2, \\
V^μ_ν(t) = V^μ_{ν0} + δ^μ_νT_2 t, \hspace 1em U_ν = U^0_ν + U^1_ν t,
$$
where the coefficients $X^μ_j$, $T_i$, $V^μ_{ν0}$, $U^i_ν$ for $i = 0, 1$ and $j = 0, 1, 2$ are all constants.
From this, we get the following solution for the field $δ$:
$$
δt = T_0 + T_1 t + T_2 t^2 + U^0_νx^ν + U^1_νx^νt, \\
δx^μ = X^μ_0 + X^μ_1 t + V^μ_{ν0}x^ν + T_2tx^μ + U^1_νx^νx^μ.
$$
Characterization of the Symmetry Group
Applying the condition $x^0 = t$, we obtain the following identities:
$$X^0_0 = T_0, \hspace 1em X^0_1 = T_1, \hspace 1em V^0_{ν0} = U^0_ν, \hspace 1em T_2= U^1_0,$$
and we can eliminate the transforms for $δt$ entirely, since they are already included in those for $δx^0$. Thus we find the following:
$$δ = \left(U^1_ν + δ^0_νT_2\right)\left(x^νx^ρ\frac{∂}{∂x^ρ}\right) + X^μ_0\frac{∂}{∂x^μ} + \left(V^μ_{ν0} + δ^0_νX^μ_1\right)\left(x^ν\frac{∂}{∂x^μ}\right).$$
Defining the dilation operator
$$D ≡ x^μ\frac{∂}{∂x^μ}$$
the independent vector fields comprising the symmetries may be written as:
$$P_μ = \frac{∂}{∂x^μ}, \hspace 1em Q^ν = x^νD, \hspace 1em R^ν_ρ = x^ν\frac{∂}{∂x^ρ},$$
and the infinitesimal transformations as
$$δx^μ = a^μ + b_νx^νx^μ + c^μ_νx^ν.$$
As can be seen by the following bracket relations,
$$\begin{array}
\hline
\left[\_,\_\right] &│& P_ρ & Q^σ & R^σ_ρ & D \\
\hline
P_μ &│& 0 & δ^σ_μD + R^σ_μ & δ^σ_μP_ρ & P_μ \\
Q^ν &│& -δ^ν_ρD - R^ν_ρ & 0 & -δ^ν_ρQ^σ & -Q^ν \\
R^ν_μ &│& -δ^ν_ρP_μ & δ^σ_μQ^ν & δ^σ_μR^ν_ρ - δ^ν_ρR^σ_μ & 0 \\
D &│& -P_ρ & Q^σ & 0 & 0
\end{array}$$
this is just the Lie algebra for $\text{SL}(n+1)$.
The symmetry group can be generalized to $\text{GL}(n+1)$ by removing the condition $D = R^ρ_ρ$ and treating $D$ as an independent generator. Equivalently, this may be carried out by expanding to a coordinate representation in $n+1$ dimensions, which may be obtained by adding a scaling coordinate $u$ and making the following replacements:
$$P_ρ = u\frac{∂}{∂x^ρ}, \hspace 1em Q^ν = -x^ν\frac{∂}{∂u}, \hspace 1em R^ν_ρ = x^ν\frac{∂}{∂x^ρ}, \hspace 1em D = -u\frac{∂}{∂u}.$$
The reduction to the $n$ dimensional version is recovered with the following replacements:
$$u → 1, \hspace 1em \frac{∂}{∂u} → -D.$$
