I'm not sure if there is supposed to be a relationship between the t intervals and rho, or whether the fact that there are 15 rows and 15 intervals is a coincidence. I'm assuming it's a coincidence, but hopefully you will see how to proceed from the principle of the following.
I've made up some data just so it would run.
lambda = linspace(0,2,17);
rho = rand(15,121);
t = linspace(0,1,16);
Here's brute force and ignorance, so to speak, with nested loops and collecting each (scalar) integral in a 4-D array.
f = @(t,lambda,rho)exp(-lambda*t)*rho;
Q = zeros([size(rho),numel(lambda),numel(t)-1]);
for i1 = 1:size(rho,1)
for i2 = 1:size(rho,2)
for i3 = 1:numel(lambda)
for j = 1:numel(t)-1
Q(i1,i2,i3,j) = integral(@(x)f(x,lambda(i3),rho(i1,i2)),t(j),t(j)+1);
end
end
end
end
We can improve on this by using the 'ArrayValued' option of INTEGRAL, treating rho as a matrix input to the qs function rather than a scalar.
f = @(t,lambda,rho)exp(-lambda*t)*rho;
Q = zeros([size(rho),numel(lambda),numel(t)-1]);
for i3 = 1:numel(lambda)
for j = 1:numel(t)-1
Q(:,:,i3,j) = integral(@(x)f(x,lambda(i3),rho),t(j),t(j)+1,'ArrayValued',true);
end
end
But why not let lambda be an array in qs as well? We just need to be clever so that we do every combination of rho elements and lambda elements. Here I've used an outer product of vectors and reshaped it.
f = @(t,lambda,rho)reshape(rho(:)*exp(-t*lambda(:).'),[size(rho),numel(lambda)]);
Q = zeros([size(rho),numel(lambda),numel(t)-1]);
for j = 1:numel(t)-1
Q(:,:,:,j) = integral(@(x)f(x,lambda,rho),t(j),t(j)+1,'ArrayValued',true);
end
Now if each interval of t needs only one row of rho, it's a slightly different problem. But we can do that, too.
