function [t,S] = pid_dynamics(S0, Sd, tau,  Gp, Gi, Gd)


    % solving
    % S' +tau*S' = Gp* S' + Gi*(S - Sd) +Gd*S''
    
    function fvec = pid_forces(t,y)
        % y1 =S , y2 = S'
        % y' = f(t,y)
        y1= y(1);
        y2= y(2);
        fvec(2) = 1/(tau - Gd) * ( (Gp-1)*y2 + Gi*(y1-Sd) );
        fvec(1) = y2;

        fvec = fvec.';
    end

    tspan=[0,100];

    % S' for P only i.e. derivative of S
    Sp = (Gp-1)/tau * (S0 -Gp/(Gp-1)*Sd); 
    y0=[S0, Sp];

    [t,y] = ode45(@pid_forces, tspan, y0);

    S= y(:,1);

    plot(t,S);
    xlabel('Time');
    ylabel('State');


end

