How to stabilize this control system?
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 &gt; 1000) i = 1000;
    if(i &lt; 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;
}