gpu reimann sum

%full code listing


t = 0:50000/10000:50000;
omega = 0:10000/10000:10000;
K=1
wt = t'*omega;
wt = wt'; %transposed the matrix as this was originally done in the wrong order
wtgpu = gpuArray(wt);
swtgpu = sin(wtgpu);
clear wtgpu; %no longer needed
swt = gather(swtgpu);
cthw=coth(K*omega);
%%
cothm = zeros(10001,10001); %initialize a large enough matrix
for j = 1:10001
     cothm(j,j) = cthw(j);
end
cothm(1,1)=0;
integralmatrix = (omega(2)-omega(1))*swt*cothm;
% now we sum across the matrix to integrate across omega for each time step t

k = 1:10001;
coft = arrayfun(@(n)sum(integralmatrix(n,:)),k); %i don't know how this function
%works but it works fast

%%
plot(t,coft)