function zdot = zdot(t,z,w)
global A B Bw Q R S T dt nT cmax n
i = min(nT-1,floor(t/dt)+1);                    % the current time step
wi = w(i) + (w(i+1)-w(i))/dt*(t - (i-1)*dt);    % interpolate disturbance
x = z(1:n); p = z(n+1:2*n);            % state and co-state vectors
uactive = -R\(S'*x + B'*p);        % unconstrained optimal control
[usat,lambda] = sat(x,uactive,p);  % saturated controls and Lagrange multipliers
% BVP equations for states and co-states
zdot = [ A*x        + B*usat + Bw*wi
        -Q*x - A'*p - S*usat + lambda(2)*cmax*T*usat ];