Fitting data to a system of ODE

I have a system of two ODEs depending on two parameter k,h and I want to find the best k,h to fit some data. To understand if my code works I started generating some fake data from the equations. Here is what I use to fit the data (I basically followed this post ODE fitting to dataset)

func=ParametricNDSolveValue[{x'[t]==k-x[t]/(1-x[t]^2-y[t]^2),y'[t]==h-y[t]/(1-x[t]^2-y[t]^2), x[0]==0,y[0]==0},{x,y},{t,0,3},{k,h}]
model[k_, h_][i_, t_] :=   Through[func[k, h][t], List][[i]] /; And @@ NumericQ /@ {k, h, i, t};
fitted = NonlinearModelFit[ dataset ,model[k, h][i, t],{{k,2/3}, {h,0.5}}, {i, t}] ;
Show[ListPlot[{Transpose[{time,datax}],Transpose[{time,datay}]}], Plot[ {fitted[1, t], fitted[2, t]}, {t, 0,3},PlotRange->All]]
fitted["ParameterTable"]


I generated my dataset solving the ODEs above with a specific value for k and h, and taking points only at some point in time. So I have the fake dataset that is in the form {index,time,data} (index 1 =data on x, index 2 =data on y) with time going from 0 to 3 with 0.1 steps.

{{1,0.,0.},{1,0.1,0.0316679},{1,0.2,0.0600199},{1,0.3,0.0851046},{1,0.4,0.107009},{1,0.5,0.12586},{1,0.6,0.14183},{1,0.7,0.155127},{1,0.8,0.165994},{1,0.9,0.174694},{1,1.,0.181503},{1,1.1,0.186694},{1,1.2,0.19053},{1,1.3,0.193256},{1,1.4,0.19509},{1,1.5,0.196221},{1,1.6,0.196812},{1,1.7,0.196998},{1,1.8,0.196888},{1,1.9,0.19657},{1,2.,0.196113},{1,2.1,0.19557},{1,2.2,0.19498},{1,2.3,0.194374},{1,2.4,0.193772},{1,2.5,0.193189},{1,2.6,0.192636},{1,2.7,0.192118},{1,2.8,0.191639},{1,2.9,0.191199},{1,3.,0.190799},{2,0.,0.},{2,0.1,0.0475773},{2,0.2,0.0905796},{2,0.3,0.129347},{2,0.4,0.164163},{2,0.5,0.195286},{2,0.6,0.222959},{2,0.7,0.247429},{2,0.8,0.268943},{2,0.9,0.287756},{2,1.,0.304121},{2,1.1,0.318291},{2,1.2,0.330508},{2,1.3,0.341005},{2,1.4,0.349996},{2,1.5,0.35768},{2,1.6,0.364233},{2,1.7,0.369814},{2,1.8,0.374562},{2,1.9,0.378599},{2,2.,0.38203},{2,2.1,0.384946},{2,2.2,0.387423},{2,2.3,0.389528},{2,2.4,0.391318},{2,2.5,0.39284},{2,2.6,0.394135},{2,2.7,0.395238},{2,2.8,0.396177},{2,2.9,0.396977},{2,3.,0.39766}}


The code works but the fit is not completely perfected as I expected since I generated the data directly from equations (here is a picture).

i.

I used k=1/3 and h=1/2 to generate data but it gives me as output k=0.26 and h =0.5018

If you want to run the code, here is the entire code:

f=1/3;
g=0.5;
sol = NDSolve[{x'[t]==f-x[t]/(1-x[t]^2-y[t]),y'[t]==g-y[t]/(1-x[t]^2-y[t]^2),x[0]==0,y[0]==0},{x,y},{t,0,3}]
datax=Flatten[Table[x[i]/.sol,{i,0,3,0.1}]];
datay=Flatten[Table[y[i]/.sol,{i,0,3,0.1}]];
time = Range[0,3,0.1];
datasetx=Transpose[{ConstantArray[1,31],time,datax}];
datasety=Transpose[{ConstantArray[2,31],time,datay}];
dataset=Join[datasetx,datasety]
func=ParametricNDSolveValue[{x'[t]==k-x[t]/(1-x[t]^2-y[t]^2),y'[t]==h-y[t]/(1-x[t]^2-y[t]^2), x[0]==0.0,y[0]==0.0},{x,y},{t,0,3},{k,h}]
model[k_, h_][i_, t_] :=   Through[func[k, h][t], List][[i]] /; And @@ NumericQ /@ {k, h, i, t};
fitted = NonlinearModelFit[ dataset ,model[k, h][i, t],{{k,0.3}, {h,0.4}}, {i, t}] ;
Show[ListPlot[{Transpose[{time,datax}],Transpose[{time,datay}]}], Plot[ {fitted[1, t], fitted[2, t]}, {t, 0,3},PlotRange->All]]
fitted["ParameterTable"]

• Please provide your dataset too! Commented Dec 16, 2021 at 14:45
• @UlrichNeumann I just edited the question adding more details on how i created the fake dataset. Commented Dec 16, 2021 at 14:52
• I would expect dataset in the form {{t[i],{x[i],y[i]}},...}, but your code doesn't run. Commented Dec 16, 2021 at 15:00
• As I said i followed what is done here mathematica.stackexchange.com/questions/84063/… which is what I need to do. I defined the model and the two index are used to refer to equation 1 or 2. I edited the question with the full code at the end that should run now. Commented Dec 16, 2021 at 15:11

dataset = Join[Table[{1, t, func[1/3, 1/2][[1]][t]}, {t, 0, 3, .1}],