function [elxp,elqp]=solve_elasticities(t,x,p,x0,dxdp,q,dqdparam,xdim,kdim,iplot)
%
%   ***   function [elxp,elqp]=solve_elasticities(t,x,p,x0,dxdp,q,dqdparam,xdim,kdim)   ***
%
%

% Elasticities of the variables and QoI wrt parameters
nt=length(t);
elxp=zeros(xdim,kdim,nt);
elqp=zeros(kdim+xdim,nt);

for it=1:nt
    for k=1:kdim
        for i=1:xdim
            if abs(x(i,it)) > 1.0E-16
                elxp(i,k,it)=p(k)/x(i,it)*dxdp(i,k,it);
            else
                if(iplot.elasticities == 1 && (iplot.dxdp.var == 1 || iplot.dxdp.param == 1))
%                     fprintf('Elasticity of variables: Attempted division by zero, it=%5i, k=%5i, i=%5i \n', it, k, i)
                end
            end
        end
        if abs(q(it)) > 1.0E-16
            elqp(k,it)=p(k)/q(it)*dqdparam(k,it);
        else
            if(iplot.elasticities == 1 && iplot.dqdp == 1)
%                 fprintf('Elasticity of QoI: Attempted division by zero, it=%5i, k=%5i,  \n', it, k)
            end
        end
    end
    
% Elasticities of the variables and QoI wrt initial conditions
     for j=1:xdim
        for i=1:xdim
            if abs(x(i,it)) > 1.0E-16
                elxp(i,kdim+j,it)=x0(j)/x(i,it)*dxdp(i,kdim+j,it);
            else
                if(iplot.elasticities == 1 && iplot.dxdz == 1)
%                     fprintf('Elasticity of variables: Attempted division by zero, it=%5i, j=%5i, i=%5i \n', it, j, i)
                end
            end
        end
        if abs(q(it)) > 1.0E-16
            elqp(kdim+j,it)=x0(j)/q(it)*dqdparam(kdim+j,it);
        else
            if(iplot.elasticities == 1 && iplot.dqdz == 1)
%                 fprintf('Elasticity of QoI: Attempted division by zero, it=%5i, k=%5i,  \n', it, j)
            end
        end
    end   
  
end


end