Fit data to a model of two-variable differential equations

I am trying to fit some data to the solution of a differential equation. I am using "ParametricNDSolve" to create a model which should fit the data using "FindFit".

Now, the code I propose is very similar to the examples I have seen to solve similar problems. However, the equation of this case is a little more complicated that the examples of Mathematica manuals (because is a function of two variables, and the initials conditions involves two numeric integrals).

The code keeps running for hours, and it does not appear to solve the problems at all:

data = {{21/100, 0.260882276365091`}, {6/25, 0.330910580727009`}, {27/ 100, 0.271601829283246`}, {3/10, 0.294039749043066`}, {33/100,  0.373363616981994`}, {9/25, 0.467495450916`}, {39/100,  0.50972503848512`}, {21/50, 0.639300915026114`}, {9/20,  0.679672806174314`}, {12/25, 0.693446859556429`}, {51/100,  0.697043731283207`}, {27/50, 0.736147748563448`}, {57/100,  0.706580484286456`}, {3/5, 0.719128578284264`}, {63/100,  0.762276173639189`}, {33/50, 0.788753648395851`}, {69/100,  0.802107836002146`}, {18/25, 0.803992011056852`}, {3/4,  0.796349018322968`}, {39/50, 0.795845629380594`}, {81/100,  0.808805024532776`}, {21/25, 0.812210415525446`}, {87/100,  0.811677706096654`}, {9/10, 0.796614634731033`}, {93/100,  0.802573366818874`}, {24/25, 0.801851431134633`}, {99/100,  0.799762250401537`}, {51/50, 0.803234757269179`}, {21/20,  0.812346679044962`}, {27/25, 0.804974524275207`}, {111/100,  0.813054479870107`}, {57/50, 0.802494974488922`}, {117/100,  0.811124462792484`}, {6/5, 0.810259975308694`}, {123/100,  0.810910967136978`}, {63/50, 0.815891670404585`}, {33/25,  0.803843274767398`}, {69/50, 0.814778081357444`}, {36/25,  0.809686831953762`}, {3/2, 0.829532613157218`}, {9/5,  0.815351190536115`}, {21/10, 0.828143435644695`}, {12/5,  0.791269047273677`}, {27/10, 0.821321230878138`}, {3,  0.81153212440728`}, {33/10, 0.821798500231733`}, {18/5,  0.819594302125574`}, {39/10, 0.828629593518031`}};  model = ParametricNDSolveValue[{ D[v[x, t], t] == d1*D[v[x, t], x, x], v[-5, t] == 0, v[10, t] == flux, v[x, 0] == flux*NIntegrate[Cos[y]^n1, {y, ArcTan[-b1*(x + a1)], Pi/2}]/NIntegrate[Cos[y]^n1, {y, -Pi/2, Pi/2}]}, v, {x, -5, 10}, {t, 0, 2}, {d1, flux, n1, b1, a1}]  fit = FindFit[data,  (c*model[d1, flux, n1, b1, a1][x, 2])/((c - 1)*model[d1, flux, n1, b1, a1][x, 2] + 1), {{d1, 0.001}, {flux, 0.57}, {c, 3.2}, {n1, 2}, {b1,1}, {a1, -0.4}}, x] 

I guess that this code is computationally inefficient, but I need help to understand why.

Thank you all very much

Ps: english is obviously not my first lenguage, so I apologize for the grammar mistakes.