Neat :)
Are you intending to break out the ICP pins? Generally the first thing I do :) Also a resistor on MCLR?

Once upon a time I found this paper...

http://www.semifluid.com/classes/PHY...Controller.pdf

the source contains code for a simple integer based PID loop, though it only controls one way it makes for a good start if you need it :)

I modified the code to control both ways and tweaked a few things while playing with it, I have no idea what state I left the code in but it's yours if you want it :)

Code:
// Variables and constants for P.I.D.
// -----------------------------------------------------------------------------
BYTE thePower;
short Tcurrent, Tlast, Tsetpoint;
short iState;
BYTE  iStateInc;
//#define iMax 254           // ctd was 1000.00 // Maximum iState
short iMax;
BYTE  pGain, iGain, dGain;

short temp;

#define ResolutionScale 16
#define iMin 1.00 // Minimum iState
#define tempMax (60 * ResolutionScale) // Maximum temperature  ctd temp * sensor resolution. Not a problem if to large a value ?
#define OPScaleFactor 128              // shift friendly
#define MaxOP 99

void InitPid( void )
{
 Tcurrent = 0;
 Tlast = 0;
 iState = 0;

 pGain = 200; //3;     // for 12 bit resolution keep <= 16 for extra safety
 iGain = 10;     // was 5, see iStateInc too
 dGain = 1000; //1;     // for 12 bit resolution keep <= 16  for extra safety

 iMax = (OPScaleFactor * MaxOP) / iGain;    // 10000 yields max 100% op, no need for more ?

 iStateInc = 4;  // new var, was const 4.0
}


void PID( void )
{
//void calculatePID(){
 /* The calculatePID() function calculates the error by subtracting the
 desired temperature (Tsetpoint) from the current temperature (Tcurrent).
 It then uses the calculated error to determine the proportional term,
 which is dependent solely upon the error, the integral term, which is an
 accumulation of the error and is used to clean up the final output, and
 the derivative term, which is dependent upon the change in the system.
 The function then adds the proportional and integral terms then subtracts
 the derivative term to come to a power estimate.
*/
/* // ctd
*  Derived from PHYS 315 - Cholewiak, Steven - PIC18F252 P.I.D. Heater Controller.pdf
*  Error increment is capped by the use of iStateInc so immune to bad readings
*  One way control only. Changed to 2 way
*  Easy to protect from overflow
*  Nasty jump at pv = sp, ok if tuned for max output ?
*  Slow wind down due to 'iState = iState - 1.00;', changed to be symmetrical
*/
short tempPower;
short pTerm, iTerm, dTerm, theError;

    Tsetpoint = 0x136;
    Tcurrent = temp;
    
   theError = Tcurrent - Tsetpoint;   // ctd swapped these over

   pTerm = pGain * theError;
   dTerm = dGain * (Tcurrent-Tlast);  // one blip per change, hardly worth it for slow systems imho

   if( theError > 0 ) {   // ctd was >=0
        iState = iState + iStateInc;
   }
   else {
        iState = iState - iStateInc;  // ctd was 1.00;
   }

   if( iState > iMax ) iState = iMax;
   if( iState < iMin ) iState = iMin;
   iTerm = iState * iGain;
   tempPower = pTerm + iTerm - dTerm;
   tempPower = tempPower / OPScaleFactor;
   if( tempPower >= 0 ) {
    thePower = (BYTE) tempPower;
    if( thePower < 30 )
      thePower = 30;
   }
   else
    thePower = 30;

//   if( Tcurrent < Tsetpoint )
//    thePower = 30;

   if( thePower > MaxOP ) thePower = MaxOP;
//   if( Tcurrent > tempMax ) thePower=0;

   Tlast = Tcurrent;

//     serial_print_hex8( thePower );  

   
    INTCONbits.GIE = 0;     
    PWMDuty1Reset = thePower;
    INTCONbits.GIE = 1;

    
//130: printf("%f,%f,%lu,%f\n\r", Tcurrent, Tsetpoint, thePower, tempPower);
// }

}

I suspect the problem with an integer control loop may be the limited resolution, still, worth a try.