@Markiyan Hirnyk Well, this was the old code, where I just needed to solve for one value for the parameters. What I'm trying to find out is an easy way to solve for several values of theta0, beta0 and thetaR without having to run dsolve every time. Then, somehow, I could call these solutions and extract stuff from then. 

I mean like a for loop

for beta from 1 to 100 do

 for theta0 from 1 to 100 do

    for thetaR from 1 to 100 do

       dsolve(system with initial conditions)

       plot(whatever)

       calculate all values I want 

    end do etc. etc.

then I'd have a file or a worksheet with the plots, and the calculations I want. But I'm not able to do that. The best thing I found is using the parameters option on dsolve and then using each value for the parameters I want, but this takes a lot of time since I have to calculate everything every time.

Hope this made things clearer

Here is an example for some specific values of the parameters, theta0, beta0 and thetaR. I removed some unecessary parts. I need to be able to calculate values for all functions in specific points for a family of parameters, and solve the equations I left on the worksheet. Don't know if this is possible though, any help is appreciated.

Download fermions-cutoff.mw

Here is an example for some specific values of the parameters, theta0, beta0 and thetaR. I removed some unecessary parts. I need to be able to calculate values for all functions in specific points for a family of parameters, and solve the equations I left on the worksheet. Don't know if this is possible though, any help is appreciated.

Download fermions-cutoff.mw