I want to graph ##H## vs. ##φ## but there is a ##\dot φ## and I know this is a differential equation, can somebody help me what to do here?

The only way is to solve the differential equation so that you know ##\phi(t)##. For that, you'll need more than just this equation, as you have an unknown function of time on both sides of the equation. You'd need to make use of a second equation, possibly the second Friedmann equation, to resolve the discrepancy.

The only way is to solve the differential equation so that you know ##\phi(t)##. For that, you'll need more than just this equation, as you have an unknown function of time on both sides of the equation. You'd need to make use of a second equation, possibly the second Friedmann equation, to resolve the discrepancy.

I was hoping to get rid of ##\dot φ## but it seems I can't find any relationship for that. If I'm to use the second Friedmann equation,
##\frac{\ddot a}{a} = -\frac{1}{6M_p^2}(ρ+3p)~~~~~~~~~ ^*~H = \frac{\dot a}{a}~~→~~\dot H = \frac{\ddot a}{a} - (\frac{\dot a}{a})^2~~→~~\dot H = \frac{\ddot a}{a} - H^2##

##\dot H + H = -\frac{1}{6M_p^2}(ρ+3p)##

The problem is the form of ##ρ## and ##p##. For warm inflation, should I consider ##ρ = ρ_λ + ρ_r## and ##p = p_λ + p_r##? Given that ##p_λ = -ρ_λ## and ##p_r = \frac{1}{3}ρ_r##

The only way is to solve the differential equation so that you know ##\phi(t)##. For that, you'll need more than just this equation, as you have an unknown function of time on both sides of the equation. You'd need to make use of a second equation, possibly the second Friedmann equation, to resolve the discrepancy.

If I define ##t_H = H^{-1}## (Hubble time) then it would be just an ODE so I could use the typical numerical calculation in Mathematica?

Have you considered the Klein-Gordon equation?
$$ \ddot{\phi}+3H\dot{\phi}+\dfrac{dV}{d\phi}=0$$

That would be the case in typical inflationary scenario but in warm inflation KG equation would be modified to

##\ddot{\phi}+(3H + Γ)\dot{\phi}+\dfrac{dV}{d\phi}=0##

There is an extra dissipation term ##Γ##, which I would also need later, so that is also a problem.
Basically, I want to find the relationship of the tensor to scalar ratio ##r## with the ##Γ## dissipation term by the theoretical result ##r = 16ε## where ##ε## is the Hubble slow roll parameter, but from the equations I can see, ##H##, ##φ##, and ##\dot φ## are in the way since ##ε = -\frac{\dot H}{H^2}## so I think I can numerically calculate ##H## in terms of ##φ## in order to get different values of H to again numerically calculate ##r## in terms of ##H##.