How can I solve it in Matlab?

A solar collector lays on a roof that is pitched 30 degrees from the horizontal, and the roof faces 15 degrees west of south. Now we wish to find the instantaneous power throughout the day, as well as the total energy captured for the day. We will use extraterrestrial quantities for this problem. First, for a horizontal collector on Feb. 1st in San Diego, plot Go over the entire day.For the solar collector tilted 30 degrees from the horizontal, and facing 15 degrees west of south plot GoT on the same graph. Calculate Ho for the horizontal collector using the equation given in the book. Now, integrate under both curves to determine how much energy per unit area the two collectors would collect over the entire day. Make sure you get the same result as the equation for the horizontal collector (there is no equation given for the tilted collector). What is the ratio of energy collected for the tilted surface compared to the horizontal surface?

This is my matlab code:

clear clc beta=30; gamma=15; pheta=32.7157; n=32; B=(32)*(360/365); G_sc=1367; delta=(180/pi)*(0.006918-0.399912*cosd(B)+0.070257*sind(B)-0.006758*cosd(2*B)+0.000907*sind(2*B)-0.002697*cosd(3*B)+0.00148*sind(3*B)); omega_s=acosd(-tand(pheta)*tand(delta)); sol_time=omega_s/15; %On Feb. 1st sun rises at 6:45 am and sets at 5:15 pm. This determines %x-axis range in minutes %G for horizontal surface A=1; for k=-78.6365:78.6365 G(A)=G_sc*(1+0.033*cosd((360*n)/365))*cosd(pheta)*cosd(delta)*cosd(k)+sind(pheta)*sind(delta); A=A+1; end %G for tilted surface AA=1; for k1=-78.6365:78.6365 G1(AA)=G_sc*(1+0.033*cosd((360*n)/365))*(sind(delta)*sind(pheta)*cosd(beta)-sind(delta)*cosd(pheta)*sind(beta)*cosd(gamma)+cosd(delta)*cosd(pheta)*cosd(beta)*cosd(k1)+cosd(delta)*sind(pheta)*sind(beta)*cosd(gamma)*cosd(k1)+cosd(delta)*sind(beta)*sind(gamma)*sind(k1)); AA=AA+1; end x=linspace(-350,350,length(G)); plot(x,G) hold on x=linspace(-350,350,length(G1)); plot(x,G1) hold off xlabel('Solar Time (min)') ylabel('Solar Radiation (W/m^2)') title('Solar radiation on February 1st in San Diego') legend('Horizontal Surface','Tilted Surface','Location','northwest') %Daily solar radiation on a horizontal surface (using equation) H=((24*3600*G_sc)/pi)*(1+0.033*cosd((360*n)/365))*(cosd(pheta)*cosd(delta)*sind(78.6365)+((pi*78.6365)/180)*sind(pheta)*sind(delta)) % Daily solar radiation on both surfaces using graph H_h=trapz(x,G)*60 H_tilted=trapz(x,G1)*60

But the integral from graph gives different answer H than the formula.