ODE solver and PID gain

Hi, I have a system of 12 non-linear equation and was using ODE45 to solve the equations.

I get 12 graphs vs time and now I am wondering how I can apply gains (PID) to 3 of the graphs (or all of them is fine)

EX: u_dot = a*r +p*r+...

    v_dot = p*u +r*cos(theta)... 
    r_dot = u*cos(phi) sin(theta)... 
equations continue... and I am plotting time vs u(t) for example.

Since it is not statespace,transfer function, or in simulink, im totally lost.

Any suggestion is greatly appreciated. Thanks