Solve very stiff ODE

I'm trying simulate the temperature (T) of a material that is repeatedly hit by proton beam, by solving an ode of the form: dT/dt=P-k(T-T0). The problem is that the power term P term is zero when there's no beam, and during the short time when there's a beam it is a very large number. I'm defining P in my function to be solved as

if (t>n*period)&&(t<n*period+time)
      P=Iinj*(3*ppower)*th_gcm
  else
      P=0;
  end

And the output of the function be solved is

dTdt=(-2*sigma*f*Ac*epsilon*(T^4-(T0)^4)+P*Ac)/(c*V*rho);

To solve it, I've divided the solution into parts (pulse and no_pulse), due to the fact that the beam is very short and can be missed by the solver. I then combine the two solutions.

        pulse_to=n*period;
        pulse_tf=n*period+time;
        tspan=[pulse_to,pulse_tf];
        [t,T] = ode15s('foiltemp',tspan,T0);

        T0=T(end);
        nopulse_to=n*period+time;
        nopulse_tf=(n+1)*period;
        tspan=[nopulse_to,nopulse_tf];
        [t,T] = ode15s('foiltemp',tspan,T0);

I've tried solving this problem with ode15s but I always get the same error message:

Warning: Failure at t=5.000000e+00.  Unable to meet integration
tolerances without reducing the step size below the smallest value
allowed (1.421085e-14) at time t. 
> In ode15s at 731

I'm guessing that this is because of the stiffness of the problem, but are there any ways to work around this?