sol=NDSolve[k*D[f1[k],k]==F[f2[k],f3[k]],k*D[f2[k],k]==F[f1[k],f3[k]],k*D[f3[k],k]==F[f1[k],f2[k]],f1[1]==1,f2[1]==1,f3[1]==1,{f1[k],f2[k],f3[k]},{k,kmin,kmax}];

k[r]=Interpolation[Table[{r,NIntegrate[g[R],{R,0,r}]},{r,rmin,rmax,rstep}]];

Plot[{Evaluate[f1[k] /. sol]}, {k, kmin, kmax}]

Plot[h[r, Evaluate[Evaluate[f1[k] /. sol]/.k->k[r]]],k[r],Evaluate[Evaluate[f2[k] /. sol]/.k->k[r]],Evaluate[Evaluate[f3[k] /. sol]/.k->k[r]]], {r, rmin, rmax}]