function [CpT1,CpT2,dH] = intCpdT(A,B,C,D,E,F,G,H,T1,T2)

% Example for methane between the temperatures of 298-500K: 
% A=-0.7;B=108.477;C=-42.52;D=5.86;E=.6785;F=-76.84;G=158.71;H=-74.87;T1=298;T2=500;[CpT1, CpT2, dH] = intCpdT(A,B,C,D,E,F,G,H,T1,T2)

% Input: Parameters (A,B,C,D,E,F,G,H) for the Shomate equation from NIST Chemistry web book (http://webbook.nist.gov/chemistry/)
% Output: Enthalpy change calculated by: dH = int(Cp*dT) between temperatures T1 and T2 in Kelvin
% The program will also give you the value of the heat capacity at T1 and T2
% Caution: The program uses the trapazoid rule for the integration and
% assumes 50 trapezoids will sufice. This may not always be the case.

t = T1:(T2-T1)/50:T2;
shomate = A+B.*(t/1000)+C.*(t/1000).^2+D.*(t/1000).^3+E./(t/1000).^2;

CpT1 = A+B.*(T1/1000)+C.*(T1/1000).^2+D.*(T1/1000).^3+E./(T1/1000).^2;
CpT2 = A+B.*(T2/1000)+C.*(T2/1000).^2+D.*(T2/1000).^3+E./(T2/1000).^2;

dH = trapz(t,shomate);