Given data, how to FindFit for function that returns the Last of the "ValuesOnGrid" for InterpolatingFunctions returned by multiple NDSolve-s

I have tried both, NDSolve and ParametricNDSolve, to tackle the following problem without success. I have looked at 2 other SE posts (here and here) that seem similar to mine, but wasn’t able to resolve my problem using those. Could someone point out what I’m missing? I’d also appreciate any pointers about the deeper Wolfram Language concepts causing this issue.

The problem: I have a function f of variable x with c1 and c2 as parameters:

f[c1_,c2_,x_]:=c1^2 (1 - x c2) HeavisideTheta[c2 - x]

This function feeds the definition of the parametric model, involving an NDSolve:

model[c1_, c2_, k_] := NDSolve[
{g'[x] + (f[c1, c2, x]/k) Sin[k x + g[x]]^2 == 0, g[0] == 0}, 
g,
{x, 15/c2}]

The above NDSolve returns an InterpolatingFunction for explicit values of the arguments c1, c2 and k.

Now, the object I’m ultimately interested in is the function of k obtained by taking the last value of the InterpolatingFunction, for each value of k.

I have numeric data (Reals) in the form {{x1,y1},{x2,y2},....,{xn,yn}}. What I’d like to do is to FindFit for parameters {c1,c2} in the following sense:

FindFit[data, Last[g["ValuesOnGrid"] /. First@model[c1, c2, k]], {c1, c2}, k]

This, however, gives the error message "Endpoint 15.708/c2 in {r,15.708/c2} is not a real number". I have tried setting this problem up using ParametricNDSolve as well, but to no avail. I’ve attached a screen-shot of what I see.

(* note I needed to add a Piecewise, because HeavsideTheta is non-numeric at zero *)
f[c1_, c2_, x_] := c1^2 (1 - x c2) Piecewise[{{1, c2 - x > 0}}, 0]

sol = g /. ParametricNDSolve[{g'[x] + (f[c1, c2, x]/k) Sin[k x + g[x]]^2 == 0, 
  g[0] == 0}, g, {x, 15/c2}, {c1, c2, k}];

SeedRandom[1];
data = Table[{k, 2 k^2 - RandomReal[{-2, 2}]}, {k, 0.001, 3, .1}];

(* get the endpoint value *)
getsolk[c1_?NumericQ, c2_?NumericQ, k_?NumericQ] := Last[sol[c1, c2, k]["ValuesOnGrid"]]

fit = FindFit[data, {getsolk[c1, c2, k]}, {c1, c2}, k]

(* result: {c1 -> -123.735, c2 -> -72.2024} *)

model[c1_?NumericQ, c2_?NumericQ, k_?NumericQ] := NDSolveValue[
  {g'[x] + (f[c1, c2, x]/k) Sin[k x + g[x]]^2 == 0, g[0] == 0},
  g[15/c2],
  {x, 15/c2}
];
model[1, 2, 3]