Pump control system without arduino

Your solution is correct but an easier method (I think) is to use superposition, first we will suppress r ( i.e. we will ignore r ), then get the transfer function, then suppress d then get the other transfer function, then sum the two up to get the final transfer function.

When we suppress r will then get

$$ \frac{Y}{d} = \frac{P}{1 + PC} $$

and when we suppress d we will get $$ \frac{V}{r} = \frac{C}{1 + PC} $$

but Y = VP meaning

$$ \frac{Y}{r} = \frac{PC}{1 + PC} $$

We the sum up the two values of Y to get the final response of

$$ Y = \frac{P}{1 + PC}d + \frac{PC}{1 + PC}r $$

There is no transfer function from d to y because y depends on both d and r, you can't evaluate the value y without knowing both d and r values.

If you have huge transient error spikes, then you need a derivative part (i.e the D in PID) to handle these.

A simple PID system is not hard to write in code, it's the tuning that's the difficult part - you need to know roughly. The integral part will handle the DC offset, the proportional part the gain, and the derivative the reaction to change.

Here's some code I wrote for a temperature controlled etching tank I made with a PIC16F and a DS18B20 temp sensor ages ago (it used some nichrome wire driven via PWM wrapped round a couple of ceramic tiles for the heating element, so it had a large thermal capacity)
The results, once tuned, were an overshoot a of less than 2 degrees on power up, and a variation of less than a degree or so thereafter. It's very old code for a one off project thrown together quickly, so it can probably be improved a fair bit.
I found this document "PID without a PHD" quite helpful. Also, Matlab, Scilab or Octave can be used to simulate and tune PID systems.

/** Control variables **/

int DEADBAND = 0;   // Range where no change will be made

// PID gain values 
double Kp = 10;
double Ki = 0.02;
double Kd = 5;

// Accumulator variable for Integral calculation
double I_acc = 0;

// Returns out the proportional term
double get_prop(int error)
{
    double p;
    p = (Kp * error); 
    return p;
}

// Returns the integral term
double get_int(int error)
{
    double i;
    I_acc += (error);

    i = (double)(Ki * I_acc);
    if(i > 1000) i = 1000;
    if(i < 0) i = 0;
    return i;
}

// Returns the derivative term
double get_deriv(int error)
{
    double d;
    static int prev_error;
    d = (Kd * (error - prev_error));
    prev_error = error;
    return d;
}

double get_pid(int error)
{
    double result;
    double p, i, d;

    p = get_prop(error);
    i = get_int(error);
    d = get_deriv(error);
    result = p + i + d;
    return result;
}