function [q,dqdparam]=solve_qoi(t,x,p,dxdp,xdim,kdim)
%
%   ***   [q,dqdparam]=solve_qoi(t,x,p,dxdp,xdim,kdim)   ***
%
%   Here qoi is a function of time
%
%

[nx,nt]=size(x);
if nt ~= length(t) 
    fprintf('qoi:  Houston, we have a problem in x \n')
end
np=length(p);
if np ~= kdim+xdim
    fprintf('qoi:  Houston, we have a problem in p \n')
end

% Compute quantity of interest
q=zeros(1,nt);
dqdx=zeros(xdim,nt);
dqdp=zeros(kdim+xdim,nt);
for it=1:nt
    [q(1,it),dqdx(1:xdim,it),dqdp(1:kdim+xdim,it)]=qoi(x(1:xdim,it),p,xdim,kdim);
end

% Compute sensitivity of quantity of interest wrt parameters
dqdparam=zeros(kdim+xdim,nt);
for it=1:nt
    for k=1:kdim+xdim
        dqdparam(k,it)=0;
        for m=1:xdim
            dqdparam(k,it)=dqdparam(k,it) + dqdx(m,it)*dxdp(m,k,it);
        end
            dqdparam(k,it)=dqdparam(k,it) + dqdp(k,it);
    end
end

end
