Problems computing an integral

Hi,

I am currently trying to compute an integral between two boundaries. I copy the part of the code below, in any case the only part missing is the procedure with all the parameters. For some reason I keep getting the message "wrong number of parameters" but I cant figure out what I am missing.

The additional inputs besides the integral variable which is t_time are given.

The code looks like a nightmare at some part but is just large algebra. I left it there to show that it all depends at different points on t_time, the variable with respect to which I integrate.

Any help will be much appreciated.

Thanks.

proc (1)=hyy(t_time,tau,pars,rhat);
//local H,aAlfa,bBeta,du_aAlfa,du_bBeta;
//local kappa,sigma,alpha,beta,mu,theta;
//{aAlfa,bBeta,du_aAlfa,du_bBeta,kappa,sigma,alpha,beta,mu,theta}=make_fxs(tau,pars,0,t_time);

local epars,kappa,lambda,theta,phi,beta,alpha,x,invx,psi,sigma,npars,mu;
local phic,kappac,lambdac,thetac,betac,alphac,xc,sigmac,phiQ,lambdaQ,kappaQ,xq;
local a,b,eps,ystep,astep,range,t,a0,b0,bt,at,iota,betah,i,j,h;
local p,q,alfa1,alfa2,alfa3,dualfa1,dualfa2,dualfa3,b1,a1,d1,c1,s,n,m,u;
local aAlfa,bBeta,du_aAlfa,du_bBeta,H_h;

@ define parameters @

epars=exp(pars);

alpha=zeros(nf,1);
beta=zeros(nf,1); @ diagonal of betatilde @
sigma=eye(nf); @ diag elements of sigma ==1 @
mu=zeros(nf,1);

if nf==1;
kappa=epars[1];
psi=pars[3];
theta=epars[5];
npars=5;
if model==1;
alpha=epars[npars+1:npars+2]; npars=npars+2;
beta=epars[npars+1:npars+2]; npars=npars+2;
if correlat; sigma[1,2]=pars[npars+1]; sigma[2,1]=pars[npars+2]; npars=npars+2; endif;
elseif model==2;
alpha=epars[npars+2]; npars=npars+3;
if correlat; sigma[1,2]=pars[npars+1]; npars=npars+1; endif;
endif;
endif;

@ parameters for rotated factors inv(sigmac)*f @
lambdaQ=inv(sigma)*diagrv(eye(nf),-kappa)*sigma - diagrv(eye(nf),psi.*beta);
phiQ=-psi.*alpha;

@ solutions for the Riccati equations @
//n=lambdaQ; m=-phiQ

q = -1;
p = sigma;
s = t_time; //I need to change this one each time I call "make_fxs"
n=lambdaQ;
m=phiQ;
u=0;

c1 = (-n + sqrt(n.^2 - 2.*p.*q))./(2.*q);
d1 = (1 - c1.*u).*(n + p.*u + sqrt((n+p.*u).^2 - p.*(p.*u.^2 + 2.*n.*u + 2.*q)))./(2.*n.*u + p.*u.^2 + 2.*q);
a1 = (d1 + c1).*u - 1;
b1 = (d1.*(n + 2.*q.*c1) + a1.*(n.*c1 + p))./(a1.*c1 - d1);

bBeta = (1 + a1.*exp(b1.*s))./(c1 + d1.*exp(b1.*s)); //psi
alfa1 = (m.*(a1.*c1 - d1)./(b1.*c1.*d1)); //phi
alfa2 = log((c1 + d1.*exp(b1.*s))./(c1 + d1));
alfa3 = m.*s./c1;
aAlfa = alfa1*alfa2 + alfa3;

dualfa1 = 2.*m.*q.*(1./2./q.*(u.*(1./2./q.*(n - (n.^2 - 2.*p.*q).^(1./2)) - (1./2./q.*u.*(n - (n.^2 - 2.*p.*q).^(1./2)) + 1)./(p.*u.^2 + 2.*n.*u + 2.*q).*(n + ((n + p.*u).^2 - p.*(p.*u.^2 + 2.*n.*u + 2.*q)).^(1./2) + p.*u)) + 1).*(n - (n.^2 - 2.*p.*q).^(1./2)) - (1./2./q.*u.*(n - (n.^2 - 2.*p.*q).^(1./2)) + 1)./(p.*u.^2 + 2.*n.*u + 2.*q).*(n + ((n + p.*u).^2 - p.*(p.*u.^2 + 2.*n.*u + 2.*q)).^(1./2) + p.*u)).^2.*(2.*n + 2.*p.*u)./((p - 1./2.*n./q.*(n - (n.^2 - 2.*p.*q).^(1./2))).*(u.*(1./2./q.*(n - (n.^2 - 2.*p.*q).^(1./2)) - (1./2./q.*u.*(n - (n.^2 - 2.*p.*q).^(1./2)) + 1)./(p.*u.^2 + 2.*n.*u + 2.*q).*(n + ((n + p.*u).^2 - p.*(p.*u.^2 + 2.*n.*u + 2.*q)).^(1./2) + p.*u)) + 1) - (n.^2 - 2.*p.*q).^(1./2).*(1./2./q.*u.*(n - (n.^2 - 2.*p.*q).^(1./2)) + 1)./(p.*u.^2 + 2.*n.*u + 2.*q).*(n + ((n + p.*u).^2 - p.*(p.*u.^2 + 2.*n.*u + 2.*q)).^(1./2) + p.*u))./(n - (n.^2 - 2.*p.*q).^(1./2))./(1./2./q.*u.*(n - (n.^2 - 2.*p.*q).^(1./2)) + 1)./(n + ((n + p.*u).^2 - p.*(p.*u.^2 + 2.*n.*u + 2.*q)).^(1./2) + p.*u) - m.*(1./2./q.*(u.*(1./2./q.*(n - (n.^2 - 2.*p.*q).^(1./2)) - (1./2./q.*u.*(n - (n.^2 - 2.*p.*q).^(1./2)) + 1)./(p.*u.^2 + 2.*n.*u + 2.*q).*(n + ((n + p.*u).^2 - p.*(p.*u.^2 + 2.*n.*u + 2.*q)).^(1./2) + p.*u)) + 1).*(n - (n.^2 - 2.*p.*q).^(1./2)) - (1./2./q.*u.*(n - (n.^2 - 2.*p.*q).^(1./2)) + 1)./(p.*u.^2 + 2.*n.*u + 2.*q).*(n + ((n + p.*u).^2 - p.*(p.*u.^2 + 2.*n.*u + 2.*q)).^(1./2) + p.*u)).^2./((p - 1./2.*n./q.*(n - (n.^2 - 2.*p.*q).^(1./2))).*(u.*(1./2./q.*(n - (n.^2 - 2.*p.*q).^(1./2)) - (1./2./q.*u.*(n - (n.^2 - 2.*p.*q).^(1./2)) + 1)./(p.*u.^2 + 2.*n.*u + 2.*q).*(n + ((n + p.*u).^2 - p.*(p.*u.^2 + 2.*n.*u + 2.*q)).^(1./2) + p.*u)) + 1) - (n.^2 - 2.*p.*q).^(1./2).*(1./2./q.*u.*(n - (n.^2 - 2.*p.*q).^(1./2)) + 1)./(p.*u.^2 + 2.*n.*u + 2.*q).*(n + ((n + p.*u).^2 - p.*(p.*u.^2 + 2.*n.*u + 2.*q)).^(1./2) + p.*u))./(1./2./q.*u.*(n - (n.^2 - 2.*p.*q).^(1./2)) + 1).^2.*(p.*u.^2 + 2.*n.*u + 2.*q)./(n + ((n + p.*u).^2 - p.*(p.*u.^2 + 2.*n.*u + 2.*q)).^(1./2) + p.*u) + 2.*m.*q.*(1./2./q.*(u.*(1./2./q.*(n - (n.^2 - 2.*p.*q).^(1./2)) - (1./2./q.*u.*(n - (n.^2 - 2.*p.*q).^(1./2)) + 1)./(p.*u.^2 + 2.*n.*u + 2.*q).*(n + ((n + p.*u).^2 - p.*(p.*u.^2 + 2.*n.*u + 2.*q)).^(1./2) + p.*u)) + 1).*(n - (n.^2 - 2.*p.*q).^(1./2)) - (1./2./q.*u.*(n - (n.^2 - 2.*p.*q).^(1./2)) + 1)./(p.*u.^2 + 2.*n.*u + 2.*q).*(n + ((n + p.*u).^2 - p.*(p.*u.^2 + 2.*n.*u + 2.*q)).^(1./2) + p.*u)).^2./((p - 1./2.*n./q.*(n - (n.^2 - 2.*p.*q).^(1./2))).*(u.*(1./2./q.*(n - (n.^2 - 2.*p.*q).^(1./2)) - (1./2./q.*u.*(n - (n.^2 - 2.*p.*q).^(1./2)) + 1)./(p.*u.^2 + 2.*n.*u + 2.*q).*(n + ((n + p.*u).^2 - p.*(p.*u.^2 + 2.*n.*u + 2.*q)).^(1./2) + p.*u)) + 1) - (n.^2 - 2.*p.*q).^(1./2).*(1./2./q.*u.*(n - (n.^2 - 2.*p.*q).^(1./2)) + 1)./(p.*u.^2 + 2.*n.*u + 2.*q).*(n + ((n + p.*u).^2 - p.*(p.*u.^2 + 2.*n.*u + 2.*q)).^(1./2) + p.*u)).^2./(n - (n.^2 - 2.*p.*q).^(1./2))./(1./2./q.*u.*(n - (n.^2 - 2.*p.*q).^(1./2)) + 1).*(p.*u.^2 + 2.*n.*u + 2.*q)./(n + ((n + p.*u).^2 - p.*(p.*u.^2 + 2.*n.*u + 2.*q)).^(1./2) + p.*u).*((p - 1./2.*n./q.*(n - (n.^2 - 2.*p.*q).^(1./2))).*(u.*((p - 1./2.*(p.*(2.*n + 2.*p.*u) - 2.*p.*(n + p.*u))./((n + p.*u).^2 - p.*(p.*u.^2 + 2.*n.*u + 2.*q)).^(1./2)).*(1./2./q.*u.*(n - (n.^2 - 2.*p.*q).^(1./2)) + 1)./(p.*u.^2 + 2.*n.*u + 2.*q) - (2.*n + 2.*p.*u).*(1./2./q.*u.*(n - (n.^2 - 2.*p.*q).^(1./2)) + 1)./(p.*u.^2 + 2.*n.*u + 2.*q).^2.*(n + ((n + p.*u).^2 - p.*(p.*u.^2 + 2.*n.*u + 2.*q)).^(1./2) + p.*u) + 1./2./q.*(n - (n.^2 - 2.*p.*q).^(1./2))./(p.*u.^2 + 2.*n.*u + 2.*q).*(n + ((n + p.*u).^2 - p.*(p.*u.^2 + 2.*n.*u + 2.*q)).^(1./2) + p.*u)) - 1./2./q.*(n - (n.^2 - 2.*p.*q).^(1./2)) + (1./2./q.*u.*(n - (n.^2 - 2.*p.*q).^(1./2)) + 1)./(p.*u.^2 + 2.*n.*u + 2.*q).*(n + ((n + p.*u).^2 - p.*(p.*u.^2 + 2.*n.*u + 2.*q)).^(1./2) + p.*u)) + (p - 1./2.*(p.*(2.*n + 2.*p.*u) - 2.*p.*(n + p.*u))./((n + p.*u).^2 - p.*(p.*u.^2 + 2.*n.*u + 2.*q)).^(1./2)).*(n.^2 - 2.*p.*q).^(1./2).*(1./2./q.*u.*(n - (n.^2 - 2.*p.*q).^(1./2)) + 1)./(p.*u.^2 + 2.*n.*u + 2.*q) - (2.*n + 2.*p.*u).*(n.^2 - 2.*p.*q).^(1./2).*(1./2./q.*u.*(n - (n.^2 - 2.*p.*q).^(1./2)) + 1)./(p.*u.^2 + 2.*n.*u + 2.*q).^2.*(n + ((n + p.*u).^2 - p.*(p.*u.^2 + 2.*n.*u + 2.*q)).^(1./2) + p.*u) + 1./2./q.*(n.^2 - 2.*p.*q).^(1./2).*(n - (n.^2 - 2.*p.*q).^(1./2))./(p.*u.^2 + 2.*n.*u + 2.*q).*(n + ((n + p.*u).^2 - p.*(p.*u.^2 + 2.*n.*u + 2.*q)).^(1./2) + p.*u)) - 4.*m.*q.*(1./2./q.*(u.*(1./2./q.*(n - (n.^2 - 2.*p.*q).^(1./2)) - (1./2./q.*u.*(n - (n.^2 - 2.*p.*q).^(1./2)) + 1)./(p.*u.^2 + 2.*n.*u + 2.*q).*(n + ((n + p.*u).^2 - p.*(p.*u.^2 + 2.*n.*u + 2.*q)).^(1./2) + p.*u)) + 1).*(n - (n.^2 - 2.*p.*q).^(1./2)) - (1./2./q.*u.*(n - (n.^2 - 2.*p.*q).^(1./2)) + 1)./(p.*u.^2 + 2.*n.*u + 2.*q).*(n + ((n + p.*u).^2 - p.*(p.*u.^2 + 2.*n.*u + 2.*q)).^(1./2) + p.*u))./((p - 1./2.*n./q.*(n - (n.^2 - 2.*p.*q).^(1./2))).*(u.*(1./2./q.*(n - (n.^2 - 2.*p.*q).^(1./2)) - (1./2./q.*u.*(n - (n.^2 - 2.*p.*q).^(1./2)) + 1)./(p.*u.^2 + 2.*n.*u + 2.*q).*(n + ((n + p.*u).^2 - p.*(p.*u.^2 + 2.*n.*u + 2.*q)).^(1./2) + p.*u)) + 1) - (n.^2 - 2.*p.*q).^(1./2).*(1./2./q.*u.*(n - (n.^2 - 2.*p.*q).^(1./2)) + 1)./(p.*u.^2 + 2.*n.*u + 2.*q).*(n + ((n + p.*u).^2 - p.*(p.*u.^2 + 2.*n.*u + 2.*q)).^(1./2) + p.*u))./(n - (n.^2 - 2.*p.*q).^(1./2))./(1./2./q.*u.*(n - (n.^2 - 2.*p.*q).^(1./2)) + 1).*(p.*u.^2 + 2.*n.*u + 2.*q)./(n + ((n + p.*u).^2 - p.*(p.*u.^2 + 2.*n.*u + 2.*q)).^(1./2) + p.*u).*((p - 1./2.*(p.*(2.*n + 2.*p.*u) - 2.*p.*(n + p.*u))./((n + p.*u).^2 - p.*(p.*u.^2 + 2.*n.*u + 2.*q)).^(1./2)).*(1./2./q.*u.*(n - (n.^2 - 2.*p.*q).^(1./2)) + 1)./(p.*u.^2 + 2.*n.*u + 2.*q) + 1./2./q.*(n - (n.^2 - 2.*p.*q).^(1./2)).*(u.*((p - 1./2.*(p.*(2.*n + 2.*p.*u) - 2.*p.*(n + p.*u))./((n + p.*u).^2 - p.*(p.*u.^2 + 2.*n.*u + 2.*q)).^(1./2)).*(1./2./q.*u.*(n - (n.^2 - 2.*p.*q).^(1./2)) + 1)./(p.*u.^2 + 2.*n.*u + 2.*q) - (2.*n + 2.*p.*u).*(1./2./q.*u.*(n - (n.^2 - 2.*p.*q).^(1./2)) + 1)./(p.*u.^2 + 2.*n.*u + 2.*q).^2.*(n + ((n + p.*u).^2 - p.*(p.*u.^2 + 2.*n.*u + 2.*q)).^(1./2) + p.*u) + 1./2./q.*(n - (n.^2 - 2.*p.*q).^(1./2))./(p.*u.^2 + 2.*n.*u + 2.*q).*(n + ((n + p.*u).^2 - p.*(p.*u.^2 + 2.*n.*u + 2.*q)).^(1./2) + p.*u)) - 1./2./q.*(n - (n.^2 - 2.*p.*q).^(1./2)) + (1./2./q.*u.*(n - (n.^2 - 2.*p.*q).^(1./2)) + 1)./(p.*u.^2 + 2.*n.*u + 2.*q).*(n + ((n + p.*u).^2 - p.*(p.*u.^2 + 2.*n.*u + 2.*q)).^(1./2) + p.*u)) - (2.*n + 2.*p.*u).*(1./2./q.*u.*(n - (n.^2 - 2.*p.*q).^(1./2)) + 1)./(p.*u.^2 + 2.*n.*u + 2.*q).^2.*(n + ((n + p.*u).^2 - p.*(p.*u.^2 + 2.*n.*u + 2.*q)).^(1./2) + p.*u) + 1./2./q.*(n - (n.^2 - 2.*p.*q).^(1./2))./(p.*u.^2 + 2.*n.*u + 2.*q).*(n + ((n + p.*u).^2 - p.*(p.*u.^2 + 2.*n.*u + 2.*q)).^(1./2) + p.*u)) - 2.*m.*q.*(1./2./q.*(u.*(1./2./q.*(n - (n.^2 - 2.*p.*q).^(1./2)) - (1./2./q.*u.*(n - (n.^2 - 2.*p.*q).^(1./2)) + 1)./(p.*u.^2 + 2.*n.*u + 2.*q).*(n + ((n + p.*u).^2 - p.*(p.*u.^2 + 2.*n.*u + 2.*q)).^(1./2) + p.*u)) + 1).*(n - (n.^2 - 2.*p.*q).^(1./2)) - (1./2./q.*u.*(n - (n.^2 - 2.*p.*q).^(1./2)) + 1)./(p.*u.^2 + 2.*n.*u + 2.*q).*(n + ((n + p.*u).^2 - p.*(p.*u.^2 + 2.*n.*u + 2.*q)).^(1./2) + p.*u)).^2./((p - 1./2.*n./q.*(n - (n.^2 - 2.*p.*q).^(1./2))).*(u.*(1./2./q.*(n - (n.^2 - 2.*p.*q).^(1./2)) - (1./2./q.*u.*(n - (n.^2 - 2.*p.*q).^(1./2)) + 1)./(p.*u.^2 + 2.*n.*u + 2.*q).*(n + ((n + p.*u).^2 - p.*(p.*u.^2 + 2.*n.*u + 2.*q)).^(1./2) + p.*u)) + 1) - (n.^2 - 2.*p.*q).^(1./2).*(1./2./q.*u.*(n - (n.^2 - 2.*p.*q).^(1./2)) + 1)./(p.*u.^2 + 2.*n.*u + 2.*q).*(n + ((n + p.*u).^2 - p.*(p.*u.^2 + 2.*n.*u + 2.*q)).^(1./2) + p.*u)).*(p - 1./2.*(p.*(2.*n + 2.*p.*u) - 2.*p.*(n + p.*u))./((n + p.*u).^2 - p.*(p.*u.^2 + 2.*n.*u + 2.*q)).^(1./2))./(n - (n.^2 - 2.*p.*q).^(1./2))./(1./2./q.*u.*(n - (n.^2 - 2.*p.*q).^(1./2)) + 1).*(p.*u.^2 + 2.*n.*u + 2.*q)./(n + ((n + p.*u).^2 - p.*(p.*u.^2 + 2.*n.*u + 2.*q)).^(1./2) + p.*u).^2;
dualfa2 = -(1./2./q.*(n - (n.^2 - 2.*p.*q).^(1./2)) - (1./2./q.*u.*(n - (n.^2 - 2.*p.*q).^(1./2)) + 1)./(p.*u.^2 + 2.*n.*u + 2.*q).*(n + ((n + p.*u).^2 - p.*(p.*u.^2 + 2.*n.*u + 2.*q)).^(1./2) + p.*u)).*(1./(1./2./q.*(n - (n.^2 - 2.*p.*q).^(1./2)) - (1./2./q.*u.*(n - (n.^2 - 2.*p.*q).^(1./2)) + 1)./(p.*u.^2 + 2.*n.*u + 2.*q).*(n + ((n + p.*u).^2 - p.*(p.*u.^2 + 2.*n.*u + 2.*q)).^(1./2) + p.*u)).*(exp(-s./(1./2./q.*(u.*(1./2./q.*(n - (n.^2 - 2.*p.*q).^(1./2)) - (1./2./q.*u.*(n - (n.^2 - 2.*p.*q).^(1./2)) + 1)./(p.*u.^2 + 2.*n.*u + 2.*q).*(n + ((n + p.*u).^2 - p.*(p.*u.^2 + 2.*n.*u + 2.*q)).^(1./2) + p.*u)) + 1).*(n - (n.^2 - 2.*p.*q).^(1./2)) - (1./2./q.*u.*(n - (n.^2 - 2.*p.*q).^(1./2)) + 1)./(p.*u.^2 + 2.*n.*u + 2.*q).*(n + ((n + p.*u).^2 - p.*(p.*u.^2 + 2.*n.*u + 2.*q)).^(1./2) + p.*u)).*((p - 1./2.*n./q.*(n - (n.^2 - 2.*p.*q).^(1./2))).*(u.*(1./2./q.*(n - (n.^2 - 2.*p.*q).^(1./2)) - (1./2./q.*u.*(n - (n.^2 - 2.*p.*q).^(1./2)) + 1)./(p.*u.^2 + 2.*n.*u + 2.*q).*(n + ((n + p.*u).^2 - p.*(p.*u.^2 + 2.*n.*u + 2.*q)).^(1./2) + p.*u)) + 1) - (n.^2 - 2.*p.*q).^(1./2).*(1./2./q.*u.*(n - (n.^2 - 2.*p.*q).^(1./2)) + 1)./(p.*u.^2 + 2.*n.*u + 2.*q).*(n + ((n + p.*u).^2 - p.*(p.*u.^2 + 2.*n.*u + 2.*q)).^(1./2) + p.*u))).*(p - 1./2.*(p.*(2.*n + 2.*p.*u) - 2.*p.*(n + p.*u))./((n + p.*u).^2 - p.*(p.*u.^2 + 2.*n.*u + 2.*q)).^(1./2)).*(1./2./q.*u.*(n - (n.^2 - 2.*p.*q).^(1./2)) + 1)./(p.*u.^2 + 2.*n.*u + 2.*q) + exp(-s./(1./2./q.*(u.*(1./2./q.*(n - (n.^2 - 2.*p.*q).^(1./2)) - (1./2./q.*u.*(n - (n.^2 - 2.*p.*q).^(1./2)) + 1)./(p.*u.^2 + 2.*n.*u + 2.*q).*(n + ((n + p.*u).^2 - p.*(p.*u.^2 + 2.*n.*u + 2.*q)).^(1./2) + p.*u)) + 1).*(n - (n.^2 - 2.*p.*q).^(1./2)) - (1./2./q.*u.*(n - (n.^2 - 2.*p.*q).^(1./2)) + 1)./(p.*u.^2 + 2.*n.*u + 2.*q).*(n + ((n + p.*u).^2 - p.*(p.*u.^2 + 2.*n.*u + 2.*q)).^(1./2) + p.*u)).*((p - 1./2.*n./q.*(n - (n.^2 - 2.*p.*q).^(1./2))).*(u.*(1./2./q.*(n - (n.^2 - 2.*p.*q).^(1./2)) - (1./2./q.*u.*(n - (n.^2 - 2.*p.*q).^(1./2)) + 1)./(p.*u.^2 + 2.*n.*u + 2.*q).*(n + ((n + p.*u).^2 - p.*(p.*u.^2 + 2.*n.*u + 2.*q)).^(1./2) + p.*u)) + 1) - (n.^2 - 2.*p.*q).^(1./2).*(1./2./q.*u.*(n - (n.^2 - 2.*p.*q).^(1./2)) + 1)./(p.*u.^2 + 2.*n.*u + 2.*q).*(n + ((n + p.*u).^2 - p.*(p.*u.^2 + 2.*n.*u + 2.*q)).^(1./2) + p.*u))).*(s./(1./2./q.*(u.*(1./2./q.*(n - (n.^2 - 2.*p.*q).^(1./2)) - (1./2./q.*u.*(n - (n.^2 - 2.*p.*q).^(1./2)) + 1)./(p.*u.^2 + 2.*n.*u + 2.*q).*(n + ((n + p.*u).^2 - p.*(p.*u.^2 + 2.*n.*u + 2.*q)).^(1./2) + p.*u)) + 1).*(n - (n.^2 - 2.*p.*q).^(1./2)) - (1./2./q.*u.*(n - (n.^2 - 2.*p.*q).^(1./2)) + 1)./(p.*u.^2 + 2.*n.*u + 2.*q).*(n + ((n + p.*u).^2 - p.*(p.*u.^2 + 2.*n.*u + 2.*q)).^(1./2) + p.*u)).*((p - 1./2.*n./q.*(n - (n.^2 - 2.*p.*q).^(1./2))).*(u.*((p - 1./2.*(p.*(2.*n + 2.*p.*u) - 2.*p.*(n + p.*u))./((n + p.*u).^2 - p.*(p.*u.^2 + 2.*n.*u + 2.*q)).^(1./2)).*(1./2./q.*u.*(n - (n.^2 - 2.*p.*q).^(1./2)) + 1)./(p.*u.^2 + 2.*n.*u + 2.*q) - (2.*n + 2.*p.*u).*(1./2./q.*u.*(n - (n.^2 - 2.*p.*q).^(1./2)) + 1)./(p.*u.^2 + 2.*n.*u + 2.*q).^2.*(n + ((n + p.*u).^2 - p.*(p.*u.^2 + 2.*n.*u + 2.*q)).^(1./2) + p.*u) + 1./2./q.*(n - (n.^2 - 2.*p.*q).^(1./2))./(p.*u.^2 + 2.*n.*u + 2.*q).*(n + ((n + p.*u).^2 - p.*(p.*u.^2 + 2.*n.*u + 2.*q)).^(1./2) + p.*u)) - 1./2./q.*(n - (n.^2 - 2.*p.*q).^(1./2)) + (1./2./q.*u.*(n - (n.^2 - 2.*p.*q).^(1./2)) + 1)./(p.*u.^2 + 2.*n.*u + 2.*q).*(n + ((n + p.*u).^2 - p.*(p.*u.^2 + 2.*n.*u + 2.*q)).^(1./2) + p.*u)) + (p - 1./2.*(p.*(2.*n + 2.*p.*u) - 2.*p.*(n + p.*u))./((n + p.*u).^2 - p.*(p.*u.^2 + 2.*n.*u + 2.*q)).^(1./2)).*(n.^2 - 2.*p.*q).^(1./2).*(1./2./q.*u.*(n - (n.^2 - 2.*p.*q).^(1./2)) + 1)./(p.*u.^2 + 2.*n.*u + 2.*q) - (2.*n + 2.*p.*u).*(n.^2 - 2.*p.*q).^(1./2).*(1./2./q.*u.*(n - (n.^2 - 2.*p.*q).^(1./2)) + 1)./(p.*u.^2 + 2.*n.*u + 2.*q).^2.*(n + ((n + p.*u).^2 - p.*(p.*u.^2 + 2.*n.*u + 2.*q)).^(1./2) + p.*u) + 1./2./q.*(n.^2 - 2.*p.*q).^(1./2).*(n - (n.^2 - 2.*p.*q).^(1./2))./(p.*u.^2 + 2.*n.*u + 2.*q).*(n + ((n + p.*u).^2 - p.*(p.*u.^2 + 2.*n.*u + 2.*q)).^(1./2) + p.*u)) - s./(1./2./q.*(u.*(1./2./q.*(n - (n.^2 - 2.*p.*q).^(1./2)) - (1./2./q.*u.*(n - (n.^2 - 2.*p.*q).^(1./2)) + 1)./(p.*u.^2 + 2.*n.*u + 2.*q).*(n + ((n + p.*u).^2 - p.*(p.*u.^2 + 2.*n.*u + 2.*q)).^(1./2) + p.*u)) + 1).*(n - (n.^2 - 2.*p.*q).^(1./2)) - (1./2./q.*u.*(n - (n.^2 - 2.*p.*q).^(1./2)) + 1)./(p.*u.^2 + 2.*n.*u + 2.*q).*(n + ((n + p.*u).^2 - p.*(p.*u.^2 + 2.*n.*u + 2.*q)).^(1./2) + p.*u)).^2.*((p - 1./2.*n./q.*(n - (n.^2 - 2.*p.*q).^(1./2))).*(u.*(1./2./q.*(n - (n.^2 - 2.*p.*q).^(1./2)) - (1./2./q.*u.*(n - (n.^2 - 2.*p.*q).^(1./2)) + 1)./(p.*u.^2 + 2.*n.*u + 2.*q).*(n + ((n + p.*u).^2 - p.*(p.*u.^2 + 2.*n.*u + 2.*q)).^(1./2) + p.*u)) + 1) - (n.^2 - 2.*p.*q).^(1./2).*(1./2./q.*u.*(n - (n.^2 - 2.*p.*q).^(1./2)) + 1)./(p.*u.^2 + 2.*n.*u + 2.*q).*(n + ((n + p.*u).^2 - p.*(p.*u.^2 + 2.*n.*u + 2.*q)).^(1./2) + p.*u)).*((p - 1./2.*(p.*(2.*n + 2.*p.*u) - 2.*p.*(n + p.*u))./((n + p.*u).^2 - p.*(p.*u.^2 + 2.*n.*u + 2.*q)).^(1./2)).*(1./2./q.*u.*(n - (n.^2 - 2.*p.*q).^(1./2)) + 1)./(p.*u.^2 + 2.*n.*u + 2.*q) + 1./2./q.*(n - (n.^2 - 2.*p.*q).^(1./2)).*(u.*((p - 1./2.*(p.*(2.*n + 2.*p.*u) - 2.*p.*(n + p.*u))./((n + p.*u).^2 - p.*(p.*u.^2 + 2.*n.*u + 2.*q)).^(1./2)).*(1./2./q.*u.*(n - (n.^2 - 2.*p.*q).^(1./2)) + 1)./(p.*u.^2 + 2.*n.*u + 2.*q) - (2.*n + 2.*p.*u).*(1./2./q.*u.*(n - (n.^2 - 2.*p.*q).^(1./2)) + 1)./(p.*u.^2 + 2.*n.*u + 2.*q).^2.*(n + ((n + p.*u).^2 - p.*(p.*u.^2 + 2.*n.*u + 2.*q)).^(1./2) + p.*u) + 1./2./q.*(n - (n.^2 - 2.*p.*q).^(1./2))./(p.*u.^2 + 2.*n.*u + 2.*q).*(n + ((n + p.*u).^2 - p.*(p.*u.^2 + 2.*n.*u + 2.*q)).^(1./2) + p.*u)) - 1./2./q.*(n - (n.^2 - 2.*p.*q).^(1./2)) + (1./2./q.*u.*(n - (n.^2 - 2.*p.*q).^(1./2)) + 1)./(p.*u.^2 + 2.*n.*u + 2.*q).*(n + ((n + p.*u).^2 - p.*(p.*u.^2 + 2.*n.*u + 2.*q)).^(1./2) + p.*u)) - (2.*n + 2.*p.*u).*(1./2./q.*u.*(n - (n.^2 - 2.*p.*q).^(1./2)) + 1)./(p.*u.^2 + 2.*n.*u + 2.*q).^2.*(n + ((n + p.*u).^2 - p.*(p.*u.^2 + 2.*n.*u + 2.*q)).^(1./2) + p.*u) + 1./2./q.*(n - (n.^2 - 2.*p.*q).^(1./2))./(p.*u.^2 + 2.*n.*u + 2.*q).*(n + ((n + p.*u).^2 - p.*(p.*u.^2 + 2.*n.*u + 2.*q)).^(1./2) + p.*u))).*(1./2./q.*u.*(n - (n.^2 - 2.*p.*q).^(1./2)) + 1)./(p.*u.^2 + 2.*n.*u + 2.*q).*(n + ((n + p.*u).^2 - p.*(p.*u.^2 + 2.*n.*u + 2.*q)).^(1./2) + p.*u) + 1./2./q.*exp(-s./(1./2./q.*(u.*(1./2./q.*(n - (n.^2 - 2.*p.*q).^(1./2)) - (1./2./q.*u.*(n - (n.^2 - 2.*p.*q).^(1./2)) + 1)./(p.*u.^2 + 2.*n.*u + 2.*q).*(n + ((n + p.*u).^2 - p.*(p.*u.^2 + 2.*n.*u + 2.*q)).^(1./2) + p.*u)) + 1).*(n - (n.^2 - 2.*p.*q).^(1./2)) - (1./2./q.*u.*(n - (n.^2 - 2.*p.*q).^(1./2)) + 1)./(p.*u.^2 + 2.*n.*u + 2.*q).*(n + ((n + p.*u).^2 - p.*(p.*u.^2 + 2.*n.*u + 2.*q)).^(1./2) + p.*u)).*((p - 1./2.*n./q.*(n - (n.^2 - 2.*p.*q).^(1./2))).*(u.*(1./2./q.*(n - (n.^2 - 2.*p.*q).^(1./2)) - (1./2./q.*u.*(n - (n.^2 - 2.*p.*q).^(1./2)) + 1)./(p.*u.^2 + 2.*n.*u + 2.*q).*(n + ((n + p.*u).^2 - p.*(p.*u.^2 + 2.*n.*u + 2.*q)).^(1./2) + p.*u)) + 1) - (n.^2 - 2.*p.*q).^(1./2).*(1./2./q.*u.*(n - (n.^2 - 2.*p.*q).^(1./2)) + 1)./(p.*u.^2 + 2.*n.*u + 2.*q).*(n + ((n + p.*u).^2 - p.*(p.*u.^2 + 2.*n.*u + 2.*q)).^(1./2) + p.*u))).*(n - (n.^2 - 2.*p.*q).^(1./2))./(p.*u.^2 + 2.*n.*u + 2.*q).*(n + ((n + p.*u).^2 - p.*(p.*u.^2 + 2.*n.*u + 2.*q)).^(1./2) + p.*u) - exp(-s./(1./2./q.*(u.*(1./2./q.*(n - (n.^2 - 2.*p.*q).^(1./2)) - (1./2./q.*u.*(n - (n.^2 - 2.*p.*q).^(1./2)) + 1)./(p.*u.^2 + 2.*n.*u + 2.*q).*(n + ((n + p.*u).^2 - p.*(p.*u.^2 + 2.*n.*u + 2.*q)).^(1./2) + p.*u)) + 1).*(n - (n.^2 - 2.*p.*q).^(1./2)) - (1./2./q.*u.*(n - (n.^2 - 2.*p.*q).^(1./2)) + 1)./(p.*u.^2 + 2.*n.*u + 2.*q).*(n + ((n + p.*u).^2 - p.*(p.*u.^2 + 2.*n.*u + 2.*q)).^(1./2) + p.*u)).*((p - 1./2.*n./q.*(n - (n.^2 - 2.*p.*q).^(1./2))).*(u.*(1./2./q.*(n - (n.^2 - 2.*p.*q).^(1./2)) - (1./2./q.*u.*(n - (n.^2 - 2.*p.*q).^(1./2)) + 1)./(p.*u.^2 + 2.*n.*u + 2.*q).*(n + ((n + p.*u).^2 - p.*(p.*u.^2 + 2.*n.*u + 2.*q)).^(1./2) + p.*u)) + 1) - (n.^2 - 2.*p.*q).^(1./2).*(1./2./q.*u.*(n - (n.^2 - 2.*p.*q).^(1./2)) + 1)./(p.*u.^2 + 2.*n.*u + 2.*q).*(n + ((n + p.*u).^2 - p.*(p.*u.^2 + 2.*n.*u + 2.*q)).^(1./2) + p.*u))).*(2.*n + 2.*p.*u).*(1./2./q.*u.*(n - (n.^2 - 2.*p.*q).^(1./2)) + 1)./(p.*u.^2 + 2.*n.*u + 2.*q).^2.*(n + ((n + p.*u).^2 - p.*(p.*u.^2 + 2.*n.*u + 2.*q)).^(1./2) + p.*u)) - 1./(1./2./q.*(n - (n.^2 - 2.*p.*q).^(1./2)) - (1./2./q.*u.*(n - (n.^2 - 2.*p.*q).^(1./2)) + 1)./(p.*u.^2 + 2.*n.*u + 2.*q).*(n + ((n + p.*u).^2 - p.*(p.*u.^2 + 2.*n.*u + 2.*q)).^(1./2) + p.*u)).^2.*(1./2./q.*(n - (n.^2 - 2.*p.*q).^(1./2)) - exp(-s./(1./2./q.*(u.*(1./2./q.*(n - (n.^2 - 2.*p.*q).^(1./2)) - (1./2./q.*u.*(n - (n.^2 - 2.*p.*q).^(1./2)) + 1)./(p.*u.^2 + 2.*n.*u + 2.*q).*(n + ((n + p.*u).^2 - p.*(p.*u.^2 + 2.*n.*u + 2.*q)).^(1./2) + p.*u)) + 1).*(n - (n.^2 - 2.*p.*q).^(1./2)) - (1./2./q.*u.*(n - (n.^2 - 2.*p.*q).^(1./2)) + 1)./(p.*u.^2 + 2.*n.*u + 2.*q).*(n + ((n + p.*u).^2 - p.*(p.*u.^2 + 2.*n.*u + 2.*q)).^(1./2) + p.*u)).*((p - 1./2.*n./q.*(n - (n.^2 - 2.*p.*q).^(1./2))).*(u.*(1./2./q.*(n - (n.^2 - 2.*p.*q).^(1./2)) - (1./2./q.*u.*(n - (n.^2 - 2.*p.*q).^(1./2)) + 1)./(p.*u.^2 + 2.*n.*u + 2.*q).*(n + ((n + p.*u).^2 - p.*(p.*u.^2 + 2.*n.*u + 2.*q)).^(1./2) + p.*u)) + 1) - (n.^2 - 2.*p.*q).^(1./2).*(1./2./q.*u.*(n - (n.^2 - 2.*p.*q).^(1./2)) + 1)./(p.*u.^2 + 2.*n.*u + 2.*q).*(n + ((n + p.*u).^2 - p.*(p.*u.^2 + 2.*n.*u + 2.*q)).^(1./2) + p.*u))).*(1./2./q.*u.*(n - (n.^2 - 2.*p.*q).^(1./2)) + 1)./(p.*u.^2 + 2.*n.*u + 2.*q).*(n + ((n + p.*u).^2 - p.*(p.*u.^2 + 2.*n.*u + 2.*q)).^(1./2) + p.*u)).*((p - 1./2.*(p.*(2.*n + 2.*p.*u) - 2.*p.*(n + p.*u))./((n + p.*u).^2 - p.*(p.*u.^2 + 2.*n.*u + 2.*q)).^(1./2)).*(1./2./q.*u.*(n - (n.^2 - 2.*p.*q).^(1./2)) + 1)./(p.*u.^2 + 2.*n.*u + 2.*q) - (2.*n + 2.*p.*u).*(1./2./q.*u.*(n - (n.^2 - 2.*p.*q).^(1./2)) + 1)./(p.*u.^2 + 2.*n.*u + 2.*q).^2.*(n + ((n + p.*u).^2 - p.*(p.*u.^2 + 2.*n.*u + 2.*q)).^(1./2) + p.*u) + 1./2./q.*(n - (n.^2 - 2.*p.*q).^(1./2))./(p.*u.^2 + 2.*n.*u + 2.*q).*(n + ((n + p.*u).^2 - p.*(p.*u.^2 + 2.*n.*u + 2.*q)).^(1./2) + p.*u)))./(1./2./q.*(n - (n.^2 - 2.*p.*q).^(1./2)) - exp(-s./(1./2./q.*(u.*(1./2./q.*(n - (n.^2 - 2.*p.*q).^(1./2)) - (1./2./q.*u.*(n - (n.^2 - 2.*p.*q).^(1./2)) + 1)./(p.*u.^2 + 2.*n.*u + 2.*q).*(n + ((n + p.*u).^2 - p.*(p.*u.^2 + 2.*n.*u + 2.*q)).^(1./2) + p.*u)) + 1).*(n - (n.^2 - 2.*p.*q).^(1./2)) - (1./2./q.*u.*(n - (n.^2 - 2.*p.*q).^(1./2)) + 1)./(p.*u.^2 + 2.*n.*u + 2.*q).*(n + ((n + p.*u).^2 - p.*(p.*u.^2 + 2.*n.*u + 2.*q)).^(1./2) + p.*u)).*((p - 1./2.*n./q.*(n - (n.^2 - 2.*p.*q).^(1./2))).*(u.*(1./2./q.*(n - (n.^2 - 2.*p.*q).^(1./2)) - (1./2./q.*u.*(n - (n.^2 - 2.*p.*q).^(1./2)) + 1)./(p.*u.^2 + 2.*n.*u + 2.*q).*(n + ((n + p.*u).^2 - p.*(p.*u.^2 + 2.*n.*u + 2.*q)).^(1./2) + p.*u)) + 1) - (n.^2 - 2.*p.*q).^(1./2).*(1./2./q.*u.*(n - (n.^2 - 2.*p.*q).^(1./2)) + 1)./(p.*u.^2 + 2.*n.*u + 2.*q).*(n + ((n + p.*u).^2 - p.*(p.*u.^2 + 2.*n.*u + 2.*q)).^(1./2) + p.*u))).*(1./2./q.*u.*(n - (n.^2 - 2.*p.*q).^(1./2)) + 1)./(p.*u.^2 + 2.*n.*u + 2.*q).*(n + ((n + p.*u).^2 - p.*(p.*u.^2 + 2.*n.*u + 2.*q)).^(1./2) + p.*u));
dualfa3 = 0;
du_aAlfa = dualfa1.*alfa2 + alfa1.*dualfa2 + dualfa3;

du_bBeta = (exp(-s./(1./2./q.*(u.*(1./2./q.*(n - (n.^2 - 2.*p.*q).^(1./2)) - (1./2./q.*u.*(n - (n.^2 - 2.*p.*q).^(1./2)) + 1)./(p.*u.^2 + 2.*n.*u + 2.*q).*(n + ((n + p.*u).^2 - p.*(p.*u.^2 + 2.*n.*u + 2.*q)).^(1./2) + p.*u)) + 1).*(n - (n.^2 - 2.*p.*q).^(1./2)) - (1./2./q.*u.*(n - (n.^2 - 2.*p.*q).^(1./2)) + 1)./(p.*u.^2 + 2.*n.*u + 2.*q).*(n + ((n + p.*u).^2 - p.*(p.*u.^2 + 2.*n.*u + 2.*q)).^(1./2) + p.*u)).*((p - 1./2.*n./q.*(n - (n.^2 - 2.*p.*q).^(1./2))).*(u.*(1./2./q.*(n - (n.^2 - 2.*p.*q).^(1./2)) - (1./2./q.*u.*(n - (n.^2 - 2.*p.*q).^(1./2)) + 1)./(p.*u.^2 + 2.*n.*u + 2.*q).*(n + ((n + p.*u).^2 - p.*(p.*u.^2 + 2.*n.*u + 2.*q)).^(1./2) + p.*u)) + 1) - (n.^2 - 2.*p.*q).^(1./2).*(1./2./q.*u.*(n - (n.^2 - 2.*p.*q).^(1./2)) + 1)./(p.*u.^2 + 2.*n.*u + 2.*q).*(n + ((n + p.*u).^2 - p.*(p.*u.^2 + 2.*n.*u + 2.*q)).^(1./2) + p.*u))).*(u.*(1./2./q.*(n - (n.^2 - 2.*p.*q).^(1./2)) - (1./2./q.*u.*(n - (n.^2 - 2.*p.*q).^(1./2)) + 1)./(p.*u.^2 + 2.*n.*u + 2.*q).*(n + ((n + p.*u).^2 - p.*(p.*u.^2 + 2.*n.*u + 2.*q)).^(1./2) + p.*u)) + 1) - 1)./(1./2./q.*(n - (n.^2 - 2.*p.*q).^(1./2)) - exp(-s./(1./2./q.*(u.*(1./2./q.*(n - (n.^2 - 2.*p.*q).^(1./2)) - (1./2./q.*u.*(n - (n.^2 - 2.*p.*q).^(1./2)) + 1)./(p.*u.^2 + 2.*n.*u + 2.*q).*(n + ((n + p.*u).^2 - p.*(p.*u.^2 + 2.*n.*u + 2.*q)).^(1./2) + p.*u)) + 1).*(n - (n.^2 - 2.*p.*q).^(1./2)) - (1./2./q.*u.*(n - (n.^2 - 2.*p.*q).^(1./2)) + 1)./(p.*u.^2 + 2.*n.*u + 2.*q).*(n + ((n + p.*u).^2 - p.*(p.*u.^2 + 2.*n.*u + 2.*q)).^(1./2) + p.*u)).*((p - 1./2.*n./q.*(n - (n.^2 - 2.*p.*q).^(1./2))).*(u.*(1./2./q.*(n - (n.^2 - 2.*p.*q).^(1./2)) - (1./2./q.*u.*(n - (n.^2 - 2.*p.*q).^(1./2)) + 1)./(p.*u.^2 + 2.*n.*u + 2.*q).*(n + ((n + p.*u).^2 - p.*(p.*u.^2 + 2.*n.*u + 2.*q)).^(1./2) + p.*u)) + 1) - (n.^2 - 2.*p.*q).^(1./2).*(1./2./q.*u.*(n - (n.^2 - 2.*p.*q).^(1./2)) + 1)./(p.*u.^2 + 2.*n.*u + 2.*q).*(n + ((n + p.*u).^2 - p.*(p.*u.^2 + 2.*n.*u + 2.*q)).^(1./2) + p.*u))).*(1./2./q.*u.*(n - (n.^2 - 2.*p.*q).^(1./2)) + 1)./(p.*u.^2 + 2.*n.*u + 2.*q).*(n + ((n + p.*u).^2 - p.*(p.*u.^2 + 2.*n.*u + 2.*q)).^(1./2) + p.*u)).^2.*(exp(-s./(1./2./q.*(u.*(1./2./q.*(n - (n.^2 - 2.*p.*q).^(1./2)) - (1./2./q.*u.*(n - (n.^2 - 2.*p.*q).^(1./2)) + 1)./(p.*u.^2 + 2.*n.*u + 2.*q).*(n + ((n + p.*u).^2 - p.*(p.*u.^2 + 2.*n.*u + 2.*q)).^(1./2) + p.*u)) + 1).*(n - (n.^2 - 2.*p.*q).^(1./2)) - (1./2./q.*u.*(n - (n.^2 - 2.*p.*q).^(1./2)) + 1)./(p.*u.^2 + 2.*n.*u + 2.*q).*(n + ((n + p.*u).^2 - p.*(p.*u.^2 + 2.*n.*u + 2.*q)).^(1./2) + p.*u)).*((p - 1./2.*n./q.*(n - (n.^2 - 2.*p.*q).^(1./2))).*(u.*(1./2./q.*(n - (n.^2 - 2.*p.*q).^(1./2)) - (1./2./q.*u.*(n - (n.^2 - 2.*p.*q).^(1./2)) + 1)./(p.*u.^2 + 2.*n.*u + 2.*q).*(n + ((n + p.*u).^2 - p.*(p.*u.^2 + 2.*n.*u + 2.*q)).^(1./2) + p.*u)) + 1) - (n.^2 - 2.*p.*q).^(1./2).*(1./2./q.*u.*(n - (n.^2 - 2.*p.*q).^(1./2)) + 1)./(p.*u.^2 + 2.*n.*u + 2.*q).*(n + ((n + p.*u).^2 - p.*(p.*u.^2 + 2.*n.*u + 2.*q)).^(1./2) + p.*u))).*(p - 1./2.*(p.*(2.*n + 2.*p.*u) - 2.*p.*(n + p.*u))./((n + p.*u).^2 - p.*(p.*u.^2 + 2.*n.*u + 2.*q)).^(1./2)).*(1./2./q.*u.*(n - (n.^2 - 2.*p.*q).^(1./2)) + 1)./(p.*u.^2 + 2.*n.*u + 2.*q) + exp(-s./(1./2./q.*(u.*(1./2./q.*(n - (n.^2 - 2.*p.*q).^(1./2)) - (1./2./q.*u.*(n - (n.^2 - 2.*p.*q).^(1./2)) + 1)./(p.*u.^2 + 2.*n.*u + 2.*q).*(n + ((n + p.*u).^2 - p.*(p.*u.^2 + 2.*n.*u + 2.*q)).^(1./2) + p.*u)) + 1).*(n - (n.^2 - 2.*p.*q).^(1./2)) - (1./2./q.*u.*(n - (n.^2 - 2.*p.*q).^(1./2)) + 1)./(p.*u.^2 + 2.*n.*u + 2.*q).*(n + ((n + p.*u).^2 - p.*(p.*u.^2 + 2.*n.*u + 2.*q)).^(1./2) + p.*u)).*((p - 1./2.*n./q.*(n - (n.^2 - 2.*p.*q).^(1./2))).*(u.*(1./2./q.*(n - (n.^2 - 2.*p.*q).^(1./2)) - (1./2./q.*u.*(n - (n.^2 - 2.*p.*q).^(1./2)) + 1)./(p.*u.^2 + 2.*n.*u + 2.*q).*(n + ((n + p.*u).^2 - p.*(p.*u.^2 + 2.*n.*u + 2.*q)).^(1./2) + p.*u)) + 1) - (n.^2 - 2.*p.*q).^(1./2).*(1./2./q.*u.*(n - (n.^2 - 2.*p.*q).^(1./2)) + 1)./(p.*u.^2 + 2.*n.*u + 2.*q).*(n + ((n + p.*u).^2 - p.*(p.*u.^2 + 2.*n.*u + 2.*q)).^(1./2) + p.*u))).*(s./(1./2./q.*(u.*(1./2./q.*(n - (n.^2 - 2.*p.*q).^(1./2)) - (1./2./q.*u.*(n - (n.^2 - 2.*p.*q).^(1./2)) + 1)./(p.*u.^2 + 2.*n.*u + 2.*q).*(n + ((n + p.*u).^2 - p.*(p.*u.^2 + 2.*n.*u + 2.*q)).^(1./2) + p.*u)) + 1).*(n - (n.^2 - 2.*p.*q).^(1./2)) - (1./2./q.*u.*(n - (n.^2 - 2.*p.*q).^(1./2)) + 1)./(p.*u.^2 + 2.*n.*u + 2.*q).*(n + ((n + p.*u).^2 - p.*(p.*u.^2 + 2.*n.*u + 2.*q)).^(1./2) + p.*u)).*((p - 1./2.*n./q.*(n - (n.^2 - 2.*p.*q).^(1./2))).*(u.*((p - 1./2.*(p.*(2.*n + 2.*p.*u) - 2.*p.*(n + p.*u))./((n + p.*u).^2 - p.*(p.*u.^2 + 2.*n.*u + 2.*q)).^(1./2)).*(1./2./q.*u.*(n - (n.^2 - 2.*p.*q).^(1./2)) + 1)./(p.*u.^2 + 2.*n.*u + 2.*q) - (2.*n + 2.*p.*u).*(1./2./q.*u.*(n - (n.^2 - 2.*p.*q).^(1./2)) + 1)./(p.*u.^2 + 2.*n.*u + 2.*q).^2.*(n + ((n + p.*u).^2 - p.*(p.*u.^2 + 2.*n.*u + 2.*q)).^(1./2) + p.*u) + 1./2./q.*(n - (n.^2 - 2.*p.*q).^(1./2))./(p.*u.^2 + 2.*n.*u + 2.*q).*(n + ((n + p.*u).^2 - p.*(p.*u.^2 + 2.*n.*u + 2.*q)).^(1./2) + p.*u)) - 1./2./q.*(n - (n.^2 - 2.*p.*q).^(1./2)) + (1./2./q.*u.*(n - (n.^2 - 2.*p.*q).^(1./2)) + 1)./(p.*u.^2 + 2.*n.*u + 2.*q).*(n + ((n + p.*u).^2 - p.*(p.*u.^2 + 2.*n.*u + 2.*q)).^(1./2) + p.*u)) + (p - 1./2.*(p.*(2.*n + 2.*p.*u) - 2.*p.*(n + p.*u))./((n + p.*u).^2 - p.*(p.*u.^2 + 2.*n.*u + 2.*q)).^(1./2)).*(n.^2 - 2.*p.*q).^(1./2).*(1./2./q.*u.*(n - (n.^2 - 2.*p.*q).^(1./2)) + 1)./(p.*u.^2 + 2.*n.*u + 2.*q) - (2.*n + 2.*p.*u).*(n.^2 - 2.*p.*q).^(1./2).*(1./2./q.*u.*(n - (n.^2 - 2.*p.*q).^(1./2)) + 1)./(p.*u.^2 + 2.*n.*u + 2.*q).^2.*(n + ((n + p.*u).^2 - p.*(p.*u.^2 + 2.*n.*u + 2.*q)).^(1./2) + p.*u) + 1./2./q.*(n.^2 - 2.*p.*q).^(1./2).*(n - (n.^2 - 2.*p.*q).^(1./2))./(p.*u.^2 + 2.*n.*u + 2.*q).*(n + ((n + p.*u).^2 - p.*(p.*u.^2 + 2.*n.*u + 2.*q)).^(1./2) + p.*u)) - s./(1./2./q.*(u.*(1./2./q.*(n - (n.^2 - 2.*p.*q).^(1./2)) - (1./2./q.*u.*(n - (n.^2 - 2.*p.*q).^(1./2)) + 1)./(p.*u.^2 + 2.*n.*u + 2.*q).*(n + ((n + p.*u).^2 - p.*(p.*u.^2 + 2.*n.*u + 2.*q)).^(1./2) + p.*u)) + 1).*(n - (n.^2 - 2.*p.*q).^(1./2)) - (1./2./q.*u.*(n - (n.^2 - 2.*p.*q).^(1./2)) + 1)./(p.*u.^2 + 2.*n.*u + 2.*q).*(n + ((n + p.*u).^2 - p.*(p.*u.^2 + 2.*n.*u + 2.*q)).^(1./2) + p.*u)).^2.*((p - 1./2.*n./q.*(n - (n.^2 - 2.*p.*q).^(1./2))).*(u.*(1./2./q.*(n - (n.^2 - 2.*p.*q).^(1./2)) - (1./2./q.*u.*(n - (n.^2 - 2.*p.*q).^(1./2)) + 1)./(p.*u.^2 + 2.*n.*u + 2.*q).*(n + ((n + p.*u).^2 - p.*(p.*u.^2 + 2.*n.*u + 2.*q)).^(1./2) + p.*u)) + 1) - (n.^2 - 2.*p.*q).^(1./2).*(1./2./q.*u.*(n - (n.^2 - 2.*p.*q).^(1./2)) + 1)./(p.*u.^2 + 2.*n.*u + 2.*q).*(n + ((n + p.*u).^2 - p.*(p.*u.^2 + 2.*n.*u + 2.*q)).^(1./2) + p.*u)).*((p - 1./2.*(p.*(2.*n + 2.*p.*u) - 2.*p.*(n + p.*u))./((n + p.*u).^2 - p.*(p.*u.^2 + 2.*n.*u + 2.*q)).^(1./2)).*(1./2./q.*u.*(n - (n.^2 - 2.*p.*q).^(1./2)) + 1)./(p.*u.^2 + 2.*n.*u + 2.*q) + 1./2./q.*(n - (n.^2 - 2.*p.*q).^(1./2)).*(u.*((p - 1./2.*(p.*(2.*n + 2.*p.*u) - 2.*p.*(n + p.*u))./((n + p.*u).^2 - p.*(p.*u.^2 + 2.*n.*u + 2.*q)).^(1./2)).*(1./2./q.*u.*(n - (n.^2 - 2.*p.*q).^(1./2)) + 1)./(p.*u.^2 + 2.*n.*u + 2.*q) - (2.*n + 2.*p.*u).*(1./2./q.*u.*(n - (n.^2 - 2.*p.*q).^(1./2)) + 1)./(p.*u.^2 + 2.*n.*u + 2.*q).^2.*(n + ((n + p.*u).^2 - p.*(p.*u.^2 + 2.*n.*u + 2.*q)).^(1./2) + p.*u) + 1./2./q.*(n - (n.^2 - 2.*p.*q).^(1./2))./(p.*u.^2 + 2.*n.*u + 2.*q).*(n + ((n + p.*u).^2 - p.*(p.*u.^2 + 2.*n.*u + 2.*q)).^(1./2) + p.*u)) - 1./2./q.*(n - (n.^2 - 2.*p.*q).^(1./2)) + (1./2./q.*u.*(n - (n.^2 - 2.*p.*q).^(1./2)) + 1)./(p.*u.^2 + 2.*n.*u + 2.*q).*(n + ((n + p.*u).^2 - p.*(p.*u.^2 + 2.*n.*u + 2.*q)).^(1./2) + p.*u)) - (2.*n + 2.*p.*u).*(1./2./q.*u.*(n - (n.^2 - 2.*p.*q).^(1./2)) + 1)./(p.*u.^2 + 2.*n.*u + 2.*q).^2.*(n + ((n + p.*u).^2 - p.*(p.*u.^2 + 2.*n.*u + 2.*q)).^(1./2) + p.*u) + 1./2./q.*(n - (n.^2 - 2.*p.*q).^(1./2))./(p.*u.^2 + 2.*n.*u + 2.*q).*(n + ((n + p.*u).^2 - p.*(p.*u.^2 + 2.*n.*u + 2.*q)).^(1./2) + p.*u))).*(1./2./q.*u.*(n - (n.^2 - 2.*p.*q).^(1./2)) + 1)./(p.*u.^2 + 2.*n.*u + 2.*q).*(n + ((n + p.*u).^2 - p.*(p.*u.^2 + 2.*n.*u + 2.*q)).^(1./2) + p.*u) + 1./2./q.*exp(-s./(1./2./q.*(u.*(1./2./q.*(n - (n.^2 - 2.*p.*q).^(1./2)) - (1./2./q.*u.*(n - (n.^2 - 2.*p.*q).^(1./2)) + 1)./(p.*u.^2 + 2.*n.*u + 2.*q).*(n + ((n + p.*u).^2 - p.*(p.*u.^2 + 2.*n.*u + 2.*q)).^(1./2) + p.*u)) + 1).*(n - (n.^2 - 2.*p.*q).^(1./2)) - (1./2./q.*u.*(n - (n.^2 - 2.*p.*q).^(1./2)) + 1)./(p.*u.^2 + 2.*n.*u + 2.*q).*(n + ((n + p.*u).^2 - p.*(p.*u.^2 + 2.*n.*u + 2.*q)).^(1./2) + p.*u)).*((p - 1./2.*n./q.*(n - (n.^2 - 2.*p.*q).^(1./2))).*(u.*(1./2./q.*(n - (n.^2 - 2.*p.*q).^(1./2)) - (1./2./q.*u.*(n - (n.^2 - 2.*p.*q).^(1./2)) + 1)./(p.*u.^2 + 2.*n.*u + 2.*q).*(n + ((n + p.*u).^2 - p.*(p.*u.^2 + 2.*n.*u + 2.*q)).^(1./2) + p.*u)) + 1) - (n.^2 - 2.*p.*q).^(1./2).*(1./2./q.*u.*(n - (n.^2 - 2.*p.*q).^(1./2)) + 1)./(p.*u.^2 + 2.*n.*u + 2.*q).*(n + ((n + p.*u).^2 - p.*(p.*u.^2 + 2.*n.*u + 2.*q)).^(1./2) + p.*u))).*(n - (n.^2 - 2.*p.*q).^(1./2))./(p.*u.^2 + 2.*n.*u + 2.*q).*(n + ((n + p.*u).^2 - p.*(p.*u.^2 + 2.*n.*u + 2.*q)).^(1./2) + p.*u) - exp(-s./(1./2./q.*(u.*(1./2./q.*(n - (n.^2 - 2.*p.*q).^(1./2)) - (1./2./q.*u.*(n - (n.^2 - 2.*p.*q).^(1./2)) + 1)./(p.*u.^2 + 2.*n.*u + 2.*q).*(n + ((n + p.*u).^2 - p.*(p.*u.^2 + 2.*n.*u + 2.*q)).^(1./2) + p.*u)) + 1).*(n - (n.^2 - 2.*p.*q).^(1./2)) - (1./2./q.*u.*(n - (n.^2 - 2.*p.*q).^(1./2)) + 1)./(p.*u.^2 + 2.*n.*u + 2.*q).*(n + ((n + p.*u).^2 - p.*(p.*u.^2 + 2.*n.*u + 2.*q)).^(1./2) + p.*u)).*((p - 1./2.*n./q.*(n - (n.^2 - 2.*p.*q).^(1./2))).*(u.*(1./2./q.*(n - (n.^2 - 2.*p.*q).^(1./2)) - (1./2./q.*u.*(n - (n.^2 - 2.*p.*q).^(1./2)) + 1)./(p.*u.^2 + 2.*n.*u + 2.*q).*(n + ((n + p.*u).^2 - p.*(p.*u.^2 + 2.*n.*u + 2.*q)).^(1./2) + p.*u)) + 1) - (n.^2 - 2.*p.*q).^(1./2).*(1./2./q.*u.*(n - (n.^2 - 2.*p.*q).^(1./2)) + 1)./(p.*u.^2 + 2.*n.*u + 2.*q).*(n + ((n + p.*u).^2 - p.*(p.*u.^2 + 2.*n.*u + 2.*q)).^(1./2) + p.*u))).*(2.*n + 2.*p.*u).*(1./2./q.*u.*(n - (n.^2 - 2.*p.*q).^(1./2)) + 1)./(p.*u.^2 + 2.*n.*u + 2.*q).^2.*(n + ((n + p.*u).^2 - p.*(p.*u.^2 + 2.*n.*u + 2.*q)).^(1./2) + p.*u)) - (exp(-s./(1./2./q.*(u.*(1./2./q.*(n - (n.^2 - 2.*p.*q).^(1./2)) - (1./2./q.*u.*(n - (n.^2 - 2.*p.*q).^(1./2)) + 1)./(p.*u.^2 + 2.*n.*u + 2.*q).*(n + ((n + p.*u).^2 - p.*(p.*u.^2 + 2.*n.*u + 2.*q)).^(1./2) + p.*u)) + 1).*(n - (n.^2 - 2.*p.*q).^(1./2)) - (1./2./q.*u.*(n - (n.^2 - 2.*p.*q).^(1./2)) + 1)./(p.*u.^2 + 2.*n.*u + 2.*q).*(n + ((n + p.*u).^2 - p.*(p.*u.^2 + 2.*n.*u + 2.*q)).^(1./2) + p.*u)).*((p - 1./2.*n./q.*(n - (n.^2 - 2.*p.*q).^(1./2))).*(u.*(1./2./q.*(n - (n.^2 - 2.*p.*q).^(1./2)) - (1./2./q.*u.*(n - (n.^2 - 2.*p.*q).^(1./2)) + 1)./(p.*u.^2 + 2.*n.*u + 2.*q).*(n + ((n + p.*u).^2 - p.*(p.*u.^2 + 2.*n.*u + 2.*q)).^(1./2) + p.*u)) + 1) - (n.^2 - 2.*p.*q).^(1./2).*(1./2./q.*u.*(n - (n.^2 - 2.*p.*q).^(1./2)) + 1)./(p.*u.^2 + 2.*n.*u + 2.*q).*(n + ((n + p.*u).^2 - p.*(p.*u.^2 + 2.*n.*u + 2.*q)).^(1./2) + p.*u))).*(u.*((p - 1./2.*(p.*(2.*n + 2.*p.*u) - 2.*p.*(n + p.*u))./((n + p.*u).^2 - p.*(p.*u.^2 + 2.*n.*u + 2.*q)).^(1./2)).*(1./2./q.*u.*(n - (n.^2 - 2.*p.*q).^(1./2)) + 1)./(p.*u.^2 + 2.*n.*u + 2.*q) - (2.*n + 2.*p.*u).*(1./2./q.*u.*(n - (n.^2 - 2.*p.*q).^(1./2)) + 1)./(p.*u.^2 + 2.*n.*u + 2.*q).^2.*(n + ((n + p.*u).^2 - p.*(p.*u.^2 + 2.*n.*u + 2.*q)).^(1./2) + p.*u) + 1./2./q.*(n - (n.^2 - 2.*p.*q).^(1./2))./(p.*u.^2 + 2.*n.*u + 2.*q).*(n + ((n + p.*u).^2 - p.*(p.*u.^2 + 2.*n.*u + 2.*q)).^(1./2) + p.*u)) - 1./2./q.*(n - (n.^2 - 2.*p.*q).^(1./2)) + (1./2./q.*u.*(n - (n.^2 - 2.*p.*q).^(1./2)) + 1)./(p.*u.^2 + 2.*n.*u + 2.*q).*(n + ((n + p.*u).^2 - p.*(p.*u.^2 + 2.*n.*u + 2.*q)).^(1./2) + p.*u)) - exp(-s./(1./2./q.*(u.*(1./2./q.*(n - (n.^2 - 2.*p.*q).^(1./2)) - (1./2./q.*u.*(n - (n.^2 - 2.*p.*q).^(1./2)) + 1)./(p.*u.^2 + 2.*n.*u + 2.*q).*(n + ((n + p.*u).^2 - p.*(p.*u.^2 + 2.*n.*u + 2.*q)).^(1./2) + p.*u)) + 1).*(n - (n.^2 - 2.*p.*q).^(1./2)) - (1./2./q.*u.*(n - (n.^2 - 2.*p.*q).^(1./2)) + 1)./(p.*u.^2 + 2.*n.*u + 2.*q).*(n + ((n + p.*u).^2 - p.*(p.*u.^2 + 2.*n.*u + 2.*q)).^(1./2) + p.*u)).*((p - 1./2.*n./q.*(n - (n.^2 - 2.*p.*q).^(1./2))).*(u.*(1./2./q.*(n - (n.^2 - 2.*p.*q).^(1./2)) - (1./2./q.*u.*(n - (n.^2 - 2.*p.*q).^(1./2)) + 1)./(p.*u.^2 + 2.*n.*u + 2.*q).*(n + ((n + p.*u).^2 - p.*(p.*u.^2 + 2.*n.*u + 2.*q)).^(1./2) + p.*u)) + 1) - (n.^2 - 2.*p.*q).^(1./2).*(1./2./q.*u.*(n - (n.^2 - 2.*p.*q).^(1./2)) + 1)./(p.*u.^2 + 2.*n.*u + 2.*q).*(n + ((n + p.*u).^2 - p.*(p.*u.^2 + 2.*n.*u + 2.*q)).^(1./2) + p.*u))).*(s./(1./2./q.*(u.*(1./2./q.*(n - (n.^2 - 2.*p.*q).^(1./2)) - (1./2./q.*u.*(n - (n.^2 - 2.*p.*q).^(1./2)) + 1)./(p.*u.^2 + 2.*n.*u + 2.*q).*(n + ((n + p.*u).^2 - p.*(p.*u.^2 + 2.*n.*u + 2.*q)).^(1./2) + p.*u)) + 1).*(n - (n.^2 - 2.*p.*q).^(1./2)) - (1./2./q.*u.*(n - (n.^2 - 2.*p.*q).^(1./2)) + 1)./(p.*u.^2 + 2.*n.*u + 2.*q).*(n + ((n + p.*u).^2 - p.*(p.*u.^2 + 2.*n.*u + 2.*q)).^(1./2) + p.*u)).*((p - 1./2.*n./q.*(n - (n.^2 - 2.*p.*q).^(1./2))).*(u.*((p - 1./2.*(p.*(2.*n + 2.*p.*u) - 2.*p.*(n + p.*u))./((n + p.*u).^2 - p.*(p.*u.^2 + 2.*n.*u + 2.*q)).^(1./2)).*(1./2./q.*u.*(n - (n.^2 - 2.*p.*q).^(1./2)) + 1)./(p.*u.^2 + 2.*n.*u + 2.*q) - (2.*n + 2.*p.*u).*(1./2./q.*u.*(n - (n.^2 - 2.*p.*q).^(1./2)) + 1)./(p.*u.^2 + 2.*n.*u + 2.*q).^2.*(n + ((n + p.*u).^2 - p.*(p.*u.^2 + 2.*n.*u + 2.*q)).^(1./2) + p.*u) + 1./2./q.*(n - (n.^2 - 2.*p.*q).^(1./2))./(p.*u.^2 + 2.*n.*u + 2.*q).*(n + ((n + p.*u).^2 - p.*(p.*u.^2 + 2.*n.*u + 2.*q)).^(1./2) + p.*u)) - 1./2./q.*(n - (n.^2 - 2.*p.*q).^(1./2)) + (1./2./q.*u.*(n - (n.^2 - 2.*p.*q).^(1./2)) + 1)./(p.*u.^2 + 2.*n.*u + 2.*q).*(n + ((n + p.*u).^2 - p.*(p.*u.^2 + 2.*n.*u + 2.*q)).^(1./2) + p.*u)) + (p - 1./2.*(p.*(2.*n + 2.*p.*u) - 2.*p.*(n + p.*u))./((n + p.*u).^2 - p.*(p.*u.^2 + 2.*n.*u + 2.*q)).^(1./2)).*(n.^2 - 2.*p.*q).^(1./2).*(1./2./q.*u.*(n - (n.^2 - 2.*p.*q).^(1./2)) + 1)./(p.*u.^2 + 2.*n.*u + 2.*q) - (2.*n + 2.*p.*u).*(n.^2 - 2.*p.*q).^(1./2).*(1./2./q.*u.*(n - (n.^2 - 2.*p.*q).^(1./2)) + 1)./(p.*u.^2 + 2.*n.*u + 2.*q).^2.*(n + ((n + p.*u).^2 - p.*(p.*u.^2 + 2.*n.*u + 2.*q)).^(1./2) + p.*u) + 1./2./q.*(n.^2 - 2.*p.*q).^(1./2).*(n - (n.^2 - 2.*p.*q).^(1./2))./(p.*u.^2 + 2.*n.*u + 2.*q).*(n + ((n + p.*u).^2 - p.*(p.*u.^2 + 2.*n.*u + 2.*q)).^(1./2) + p.*u)) - s./(1./2./q.*(u.*(1./2./q.*(n - (n.^2 - 2.*p.*q).^(1./2)) - (1./2./q.*u.*(n - (n.^2 - 2.*p.*q).^(1./2)) + 1)./(p.*u.^2 + 2.*n.*u + 2.*q).*(n + ((n + p.*u).^2 - p.*(p.*u.^2 + 2.*n.*u + 2.*q)).^(1./2) + p.*u)) + 1).*(n - (n.^2 - 2.*p.*q).^(1./2)) - (1./2./q.*u.*(n - (n.^2 - 2.*p.*q).^(1./2)) + 1)./(p.*u.^2 + 2.*n.*u + 2.*q).*(n + ((n + p.*u).^2 - p.*(p.*u.^2 + 2.*n.*u + 2.*q)).^(1./2) + p.*u)).^2.*((p - 1./2.*n./q.*(n - (n.^2 - 2.*p.*q).^(1./2))).*(u.*(1./2./q.*(n - (n.^2 - 2.*p.*q).^(1./2)) - (1./2./q.*u.*(n - (n.^2 - 2.*p.*q).^(1./2)) + 1)./(p.*u.^2 + 2.*n.*u + 2.*q).*(n + ((n + p.*u).^2 - p.*(p.*u.^2 + 2.*n.*u + 2.*q)).^(1./2) + p.*u)) + 1) - (n.^2 - 2.*p.*q).^(1./2).*(1./2./q.*u.*(n - (n.^2 - 2.*p.*q).^(1./2)) + 1)./(p.*u.^2 + 2.*n.*u + 2.*q).*(n + ((n + p.*u).^2 - p.*(p.*u.^2 + 2.*n.*u + 2.*q)).^(1./2) + p.*u)).*((p - 1./2.*(p.*(2.*n + 2.*p.*u) - 2.*p.*(n + p.*u))./((n + p.*u).^2 - p.*(p.*u.^2 + 2.*n.*u + 2.*q)).^(1./2)).*(1./2./q.*u.*(n - (n.^2 - 2.*p.*q).^(1./2)) + 1)./(p.*u.^2 + 2.*n.*u + 2.*q) + 1./2./q.*(n - (n.^2 - 2.*p.*q).^(1./2)).*(u.*((p - 1./2.*(p.*(2.*n + 2.*p.*u) - 2.*p.*(n + p.*u))./((n + p.*u).^2 - p.*(p.*u.^2 + 2.*n.*u + 2.*q)).^(1./2)).*(1./2./q.*u.*(n - (n.^2 - 2.*p.*q).^(1./2)) + 1)./(p.*u.^2 + 2.*n.*u + 2.*q) - (2.*n + 2.*p.*u).*(1./2./q.*u.*(n - (n.^2 - 2.*p.*q).^(1./2)) + 1)./(p.*u.^2 + 2.*n.*u + 2.*q).^2.*(n + ((n + p.*u).^2 - p.*(p.*u.^2 + 2.*n.*u + 2.*q)).^(1./2) + p.*u) + 1./2./q.*(n - (n.^2 - 2.*p.*q).^(1./2))./(p.*u.^2 + 2.*n.*u + 2.*q).*(n + ((n + p.*u).^2 - p.*(p.*u.^2 + 2.*n.*u + 2.*q)).^(1./2) + p.*u)) - 1./2./q.*(n - (n.^2 - 2.*p.*q).^(1./2)) + (1./2./q.*u.*(n - (n.^2 - 2.*p.*q).^(1./2)) + 1)./(p.*u.^2 + 2.*n.*u + 2.*q).*(n + ((n + p.*u).^2 - p.*(p.*u.^2 + 2.*n.*u + 2.*q)).^(1./2) + p.*u)) - (2.*n + 2.*p.*u).*(1./2./q.*u.*(n - (n.^2 - 2.*p.*q).^(1./2)) + 1)./(p.*u.^2 + 2.*n.*u + 2.*q).^2.*(n + ((n + p.*u).^2 - p.*(p.*u.^2 + 2.*n.*u + 2.*q)).^(1./2) + p.*u) + 1./2./q.*(n - (n.^2 - 2.*p.*q).^(1./2))./(p.*u.^2 + 2.*n.*u + 2.*q).*(n + ((n + p.*u).^2 - p.*(p.*u.^2 + 2.*n.*u + 2.*q)).^(1./2) + p.*u))).*(u.*(1./2./q.*(n - (n.^2 - 2.*p.*q).^(1./2)) - (1./2./q.*u.*(n - (n.^2 - 2.*p.*q).^(1./2)) + 1)./(p.*u.^2 + 2.*n.*u + 2.*q).*(n + ((n + p.*u).^2 - p.*(p.*u.^2 + 2.*n.*u + 2.*q)).^(1./2) + p.*u)) + 1))./(1./2./q.*(n - (n.^2 - 2.*p.*q).^(1./2)) - exp(-s./(1./2./q.*(u.*(1./2./q.*(n - (n.^2 - 2.*p.*q).^(1./2)) - (1./2./q.*u.*(n - (n.^2 - 2.*p.*q).^(1./2)) + 1)./(p.*u.^2 + 2.*n.*u + 2.*q).*(n + ((n + p.*u).^2 - p.*(p.*u.^2 + 2.*n.*u + 2.*q)).^(1./2) + p.*u)) + 1).*(n - (n.^2 - 2.*p.*q).^(1./2)) - (1./2./q.*u.*(n - (n.^2 - 2.*p.*q).^(1./2)) + 1)./(p.*u.^2 + 2.*n.*u + 2.*q).*(n + ((n + p.*u).^2 - p.*(p.*u.^2 + 2.*n.*u + 2.*q)).^(1./2) + p.*u)).*((p - 1./2.*n./q.*(n - (n.^2 - 2.*p.*q).^(1./2))).*(u.*(1./2./q.*(n - (n.^2 - 2.*p.*q).^(1./2)) - (1./2./q.*u.*(n - (n.^2 - 2.*p.*q).^(1./2)) + 1)./(p.*u.^2 + 2.*n.*u + 2.*q).*(n + ((n + p.*u).^2 - p.*(p.*u.^2 + 2.*n.*u + 2.*q)).^(1./2) + p.*u)) + 1) - (n.^2 - 2.*p.*q).^(1./2).*(1./2./q.*u.*(n - (n.^2 - 2.*p.*q).^(1./2)) + 1)./(p.*u.^2 + 2.*n.*u + 2.*q).*(n + ((n + p.*u).^2 - p.*(p.*u.^2 + 2.*n.*u + 2.*q)).^(1./2) + p.*u))).*(1./2./q.*u.*(n - (n.^2 - 2.*p.*q).^(1./2)) + 1)./(p.*u.^2 + 2.*n.*u + 2.*q).*(n + ((n + p.*u).^2 - p.*(p.*u.^2 + 2.*n.*u + 2.*q)).^(1./2) + p.*u));

H_h=exp(aAlfa+bBeta*rhat).*(du_bBeta*rhat+du_aAlfa);
retp(H_h);
endp;
let xl={1,0.6};
Lt_def=intsimp(&hyy,xl,1e-8);

17 Answers

0

The problem is that intsimp only passes one argument to your procedure. The extra arguments used by your function need to be global variables.

link

aptech

1,773

0

Thanks but I still dont get it. What are the extra arguments I am passing on to the intsimp function? I only added the bounds and the tolerance value.

Can you give me more details? is super urgent 🙁

Thanks in any case

link

Babi

6

0

The problem is that your function takes 4 inputs:

proc (1)=hyy(t_time,tau,pars,rhat);

When intsimp calls your function, it will only pass one argument to hyy. However, you cannot call your function hyy and only pass it one argument.

Let's look at a very simple example. This works:

proc (1) =  myproc(x);
    local out;
    out = sin(x);
  retp(out);
endp;
 
x_lims = { 1, 0 };
 
y = intsimp(&myproc,x_lims,1e-8);
print y;

However, if we change myproc to take more than one argument, it will not work. For example, the code below will return the error G0159 : Wrong number of parameters:

proc (1) =  myproc(x, a);
    local out;
    out = sin(x) .* a;
  retp(out);
endp;
 
x_lims = { 1, 0 };
 
y = intsimp(&myproc,x_lims,1e-8);
print y;

We can make it work, by removing a from the list of inputs to myproc and making it a global value. For example, this will work:

x_lims = { 1, 0 };

//Set global that function will use
a = 9;

y = intsimp(&myproc,x_lims,1e-8);
print y;

proc (1) =  myproc(x);
    local out;
    out = sin(x) .* a;
  retp(out);
endp;

Does that make sense?

link

aptech

1,773

0

Yes, it is perfectly clear.

But I am having now a more complicated example which I cant figure out and I get the same error of wrong number of paramenters.

Below I paste a simplified example of my code:

a=0.2;
b=0.9;
c=0.3;
d=1;

proc H_h(a,b,c,d);
local t_time;
aAlfa=a+b*t_time;
bBeta=c+d*t_time;
retp(aAlfa,bBeta);
endp;

proc hyy(t_time);
local H_h,aAlfa,bBeta,du_aAlfa,du_bBeta;
local a,b,c,d;
//local kappa,sigma,alpha,beta,mu,theta;

//{aAlfa,bBeta,du_aAlfa,du_bBeta,kappa,sigma,alpha,beta,mu,theta}=make_fxs(tau,pars,0);
{aAlfa,bBeta}=H_h(a,b,c,d);

// aAlfa=2;
// bBeta=0.5;
// du_aAlfa=0.5;
// du_bBeta=0.8;
H_h=exp(aAlfa+bBeta*t_time).*(du_bBeta*t_time+du_aAlfa);

retp(H_h);
endp;

let xl={1,0};
Lt_def=intsimp(&hyy,xl,1e-8);

 

Actually it says "Variable not initialized"

 

where is the error?

link

Babi

6

0

Here may be another way of stating my example, where again the error says Wrong Number of Parameters

 

a=0.2;
b=0.9;
c=0.3;
d=1;
rhat=3;
u=0;
tmin=0;
tmax=1;

proc alfabetaMG(a,b,c,d,t);
local aAlfa,bBeta;
aAlfa=a+b*t;
bBeta=c+d*t;
retp(aAlfa,bBeta);
endp;

proc (1) = integrandMG(t);
local H,hy,aAlfa,bBeta,du_aAlfa,du_bBeta;

{aAlfa,bBeta}=alfabetaMG(a,b,c,d,u,t);

hy=exp(aAlfa+bBeta*rhat);
H=hy;
retp(H);
endp;

proc (1) = HtMG(a,b,c,d,rhat,u,tmin,tmax);
local ret;
let xl={1,0};
ret=intsimp(&integrandMG,xl,1e-8);
retp(ret);
endp;

z1=HtMG(a,b,c,d,rhat,u,tmin,tmax);

 

 

link

Babi

6

0

This is a great question. It has some interesting concepts in it.

First error: Undefined symbol 'aAlfa'
When I run the code you pasted above, I get the error G0025 : Undefined symbol 'aAlfa'. This error is coming from the line inside of the procedure H_h:

aAlfa = a + b * t_time;

Inside of the procedure H_h, GAUSS will first see if it has a local variable named aAlfa. Since aAlfa is not declared as a local variable inside of H_h, GAUSS will next look to see if there is a global variable named aAlfa. In this case, there is not a global version either. You can fix this error by changing the procedure to either:

 
proc (2) =  H_h(a,b,c,d);
    local t_time, aAlfa, bBeta;
    aAlfa=a+b*t_time;
    bBeta=c+d*t_time;
    retp(aAlfa,bBeta);
endp;

or this:

proc (2) =  H_h(a,b,c,d);
    local t_time;
    retp(a+b*t_time, c+d*t_time);
endp;

Second error: Variable not initialized 't_time'
Now that we have fixed the first error, there is still another in the procedure H_h. The problem is that the variable t_time is declared to be a local variable in the procedure H_h. However, it is never assigned a value. When you declare a variable to be local to a procedure, you tell GAUSS that the variable will be defined inside the procedure. GAUSS will ignore any global variables with the same name.

Since it looks like you want to use the value of t_time that is passed into the procedure hyy, you should add that as an input to the procedure H_h like this:

proc (2) =  H_h(t_time,a,b,c,d);
    local aAlfa, bBeta;
    aAlfa=a+b*t_time;
    bBeta=c+d*t_time;
    retp(aAlfa,bBeta);
endp;

Notice that t_time is no longer in the local variables list, because it is passed in as an argument.
Why doesn't H_h use the variable t_time from the procedure hyy?
GAUSS procedures can use only variables that are local to the current procedure, global variables or variables that are passed in. If GAUSS procedures saw every reference to variables in other procedures, you could end up with strange bugs where your procedure sometimes had its values changed depending upon which other procedure called it.
Third error: H_h declared as variable and procedure
You cannot have a variable and a procedure with the same name in GAUSS. Make sure that you change the name of either the procedure H_h or the name of the local variable H_h.
Fourth error
It looks like you want to use the global variables a, b, c and d inside of the procedure H_h. However, the code above declares a-d to be local to the procedure hyy, but then never defines them. You can either move the definitions of a-b inside of hyy, like this:

proc (1) = hyy(t_time);
   local a, b, c, d;
    
   a = 0.2;
   b = 0.9;
   c = 0.3;
   d = 1;
   .
   .
   .
endp;

or just use the global a-d in the procedure H_t and not pass them in as inputs.

I hope this long post is clear enough to be helpful and not just long. 🙂

link

aptech

1,773

0

You are a genious, dont know how to thank you. Let me go over it in detail.

link

Babi

6

0

Aggghhhh honestly, I am going crazy. I go over all the details but this message of the Variable not initialized never goes away.

This is THE example I am trying to run

I have absolutely no clue what is the problem anymore.

Any help is highly welcomed :)))

rhat=3;

proc startval;
local a,b,c,d,rhat,u,pars;
a=0.2;
b=0.9;
c=0.3;
d=1;
u=0;
pars=a|b|c|d|u;
retp(pars);
endp;

//pars=startval;

proc (2)=alfabetaMG(pars,t);
local aAlfa,bBeta,a,b,c,d,rhat,u;
a=pars[1];
b=pars[2];
c=pars[3];
d=pars[4];
u=pars[5];
aAlfa=a+b*t;
bBeta=c+d*t;
retp(aAlfa,bBeta);
endp;

proc integrandMG(t);
local H,hy,aAlfa,bBeta;
local pars;

{aAlfa,bBeta}=alfabetaMG(pars,t);

hy=exp(aAlfa+bBeta*rhat);
retp(hy);
endp;

let xl={1,0};
z1=intsimp(&integrandMG,xl,1e-8);

link

Babi

6

0

I have a working version of the code at the bottom, but let's look at the issues first. Here is your integradMG procedure with some comments that might help:

proc integrandMG(t);
    local H,hy,aAlfa,bBeta;

    //On the line below 'pars' is declared as a
    //local variable. This tells GAUSS that
    //you will define 'pars' inside THIS PROCEDURE.
    //Further, it tells GAUSS to NOT look for a global
    //version of 'pars'
    local pars;
    
    //On the line, your procedure tries to use 'pars',
    //but the local variable 'pars' has not been defined
    {aAlfa,bBeta}=alfabetaMG(pars,t);
    
    hy=exp(aAlfa+bBeta*rhat);
    retp(hy);
endp;

Below is a working version of the code. Try to remember that variable declared local to a procedure must be defined in that procedure. If you want the procedure to use a global variable, you cannot declare it as a local variable. Finally, the variable must always be defined before it is used.

Also, in a case where you are passing a procedure to another procedure (optimization, integration, etc), it can be simpler to debug it by trying to just call your function at one point. Then pass it to the other function (maxlik, cml, optmum, intsimp, etc) after you have it working.

rhat=3;

proc startval;
    local a,b,c,d,rhat,u,pars;
    a=0.2;
    b=0.9;
    c=0.3;
    d=1;
    u=0;
    pars=a|b|c|d|u;
    retp(pars);
endp;

//Assigning the global 'pars' variable
pars=startval;

proc (2)=alfabetaMG(pars,t);
    local aAlfa,bBeta,a,b,c,d,rhat,u;
    a=pars[1];
    b=pars[2];
    c=pars[3];
    d=pars[4];
    u=pars[5];
    aAlfa=a+b*t;
    bBeta=c+d*t;
    retp(aAlfa,bBeta);
endp;

proc integrandMG(t);
    local H,hy,aAlfa,bBeta;
    
    //Using the global 'pars'
    {aAlfa,bBeta}=alfabetaMG(pars,t);
    
    hy=exp(aAlfa+bBeta*rhat);
    retp(hy);
endp;

xl = { 1, 0 };
z1 = intsimp(&integrandMG,xl,1e-8);

link

aptech

1,773

0

Now I do have a tiny new issue related to the same example I posted before.

I added one new thing and I dont know how to pass on inputs.

This is the extended example:

proc startval;
local a,b,c,d,rhat,u,pars;
a=0.2;
b=0.9;
c=0.3;
d=1;
u=0;
rhat=3;
pars=a|b|c|d|u|rhat;
retp(pars);
endp;

pars=startval;

proc (2)=alfabetaMG(pars,t);
local aAlfa,bBeta,a,b,c,d,rhat,u;
a=pars[1];
b=pars[2];
c=pars[3];
d=pars[4];
u=pars[5];
rhat=pars[6];
aAlfa=a+b*t;
bBeta=c+d*t;
retp(aAlfa,bBeta);
endp;

proc integrandMG(t);
local H,hy,aAlfa,bBeta,rhat;

{aAlfa,bBeta}=alfabetaMG(pars,t);

hy=exp(aAlfa+bBeta*rhat);
retp(hy);
endp;

proc (1)=kalman;
local rhat,z1;

rhat=3;

let xl={1,0};
z1=intsimp(&integrandMG,xl,1e-8);
retp(z1);
endp;

z1=kalman;

 

I get the error of Variable not initialized.

 

Any help?

Sorry for having so many questions.

Thanks

link

Babi

6

0

First problem
The problem is in this procedure:

proc (1)=kalman;
   local rhat,z1;

   rhat=3;

   let xl={1,0};
   z1=intsimp(&integrandMG,xl,1e-8);
   retp(z1);
endp;

I get the error of "Undefined Symobol : 'xl'". You just need to add xl to the local list for the procedure kalman, like this:

proc (1)=kalman;
   local xl, rhat,z1;

   rhat=3;

   xl={1,0};
   z1=intsimp(&integrandMG,xl,1e-8);
   retp(z1);
endp;

Also notice that I got rid of the let statement. It does not help or hurt the code, but I think it may add confusion.

Second problem
GAUSS procedures ONLY see variables that are:

1. Local to the procedure.
2. Global (i.e. defined at the main level, not defined in any GAUSS procedure)
3. Passed in to the procedure as inputs

It looks like you are trying to send the rhat that is local to kalman to your procedure integrandMG.

proc integrandMG(t);
    local H,hy,aAlfa,bBeta,rhat;
    
    {aAlfa,bBeta}=alfabetaMG(pars,t);
    
    hy=exp(aAlfa+bBeta*rhat);
    retp(hy);
endp;

proc (1)=kalman;
    local rhat,z1;
    
    //'rhat' below is a local variable to 'kalman'
    //it will not be visible in any other procedure
    //unless it is passed to that procedure
    //or defined at the global level
    rhat=3;
    
    let xl={1,0};
    z1=intsimp(&integrandMG,xl,1e-8);
    retp(z1);
endp;

link

aptech

1,773

0

Yes, the xl problem I fixed. But rhat HAS to be local to kalman.

I have rewritten the previous example and now it seems to run:

proc startval;
local a,b,c,d,u,pars;
a=0.2;
b=0.9;
c=0.3;
d=1;
u=0;
rhat=3;
pars=a|b|c|d|u;
retp(pars);
endp;

pars=startval;

proc (2)=alfabetaMG(pars,t);
local aAlfa,bBeta,a,b,c,d,u;
a=pars[1];
b=pars[2];
c=pars[3];
d=pars[4];
u=pars[5];
//rhat=pars[6];
aAlfa=a+b*t;
bBeta=c+d*t;
retp(aAlfa,bBeta);
endp;

proc integrandMG(t);
local H,hy,aAlfa,bBeta;

{aAlfa,bBeta}=alfabetaMG(pars,t);

hy=exp(aAlfa+bBeta*rhat);
retp(hy);
endp;

proc (1)=integral_test(rhat);
local xl,z1;

let xl={1,0};
z1=intsimp(&integrandMG,xl,1e-8);
retp(z1);
endp;

proc (1)=kalman;
local rhat,z1,aux;

rhat=3;

z1=integral_test(rhat);
aux=z1^2;

retp(aux);
endp;

aux=kalman;

link

Babi

6

0

No, it doesnt 🙁

I understand that Gauss only sees locals, globals and passed on arguments but since I have this procedure for the integral which does not allow to pass on anything, I have no clue how to include de rhat given that it is a local to the final procedure.

Any ideas on this?

Thanks in advance.

link

Babi

6

0

Honestly I cant see the way around it. Because rhat it is a local for the kalman procedure. That I cannot change and so I cannot pass the rhat from one process to another. But also I cannot put a proc inside another procedure and for the integral I need to express the integrand as a process, so I cannot put it directly into the kalman process.

I ran out of ideas....

link

Babi

6

0

I think this might work:

 

proc startval;
local a,b,c,d,u,pars;
a=0.2;
b=0.9;
c=0.3;
d=1;
u=0;
pars=a|b|c|d|u;
retp(pars);
endp;

pars=startval;

proc (2)=alfabetaMG(pars,t);
local aAlfa,bBeta,a,b,c,d,u;
a=pars[1];
b=pars[2];
c=pars[3];
d=pars[4];
u=pars[5];
aAlfa=a+b*t;
bBeta=c+d*t;
retp(aAlfa,bBeta);
endp;

proc integrandMG(t);
local hy,aAlfa,bBeta,rhat;

{aAlfa,bBeta}=alfabetaMG(pars,t);

hy=exp(aAlfa+bBeta*myrhat);
retp(hy);
endp;

proc (1)=integral_test(myrhat);
local xl,z1;

let xl={1,0};
z1=intsimp(&integrandMG,xl,1e-8);
retp(z1);
endp;

myrhat=8;

proc (1)=kalman;
local z1,aux,rhat;

rhat=3;
myrhat=rhat;
z1=integral_test(myrhat);
aux=z1^2;

retp(aux);
endp;

aux=kalman;

 

Could you tell me if there is any mistake in it?

sorry again and thanks!

link

Babi

6

0

Here, this version is more clear:

new;
cls;

//Set global 'pars' start value vector for use inside of 'integrandMG'
pars=startval();

//Global variable used by 'integrandMG'
myrhat=8;

//Call main procedure
aux=kalman();

proc (1)=kalman;
    local z1,aux,rhat;
    
    rhat=3;
    
    //Setting global 'myrhat' below
    //This global 'myrhat' will be referenced
    //in 'integrandMG'
    myrhat=rhat;
    z1=integral_test();
    aux=z1^2;
    
    retp(aux);
endp;

proc (1)=integral_test();
    local xl,z1;
    
    xl = { 1, 0 };
    z1 = intsimp(&integrandMG,xl,1e-8);
    retp(z1);
endp;

proc integrandMG(t);
    local hy,aAlfa,bBeta,rhat;
    
    {aAlfa,bBeta}=alfabetaMG(pars,t);
    
    hy=exp(aAlfa+bBeta*myrhat);
    retp(hy);
endp;

proc (2)=alfabetaMG(pars,t);
    local aAlfa,bBeta,a,b,c,d,u;
    a=pars[1];
    b=pars[2];
    c=pars[3];
    d=pars[4];
    u=pars[5];
    aAlfa=a+b*t;
    bBeta=c+d*t;
    retp(aAlfa,bBeta);
endp;

proc startval();
    local a,b,c,d,u,pars;
    a=0.2;
    b=0.9;
    c=0.3;
    d=1;
    u=0;
    pars=a|b|c|d|u;
    retp(pars);
endp;

And this might be even simpler, but maybe the code above is better if you need to use the procedures in another context:

new;
cls;

//Set global 'pars' start value vector for use inside of 'integrandMG'
pars=startval();

//Global variable used by 'integrandMG'
//The _g ending indicates it is a global
myrhat_g=8;

//Call main procedure
aux=kalman();

proc (1)=kalman;
    local z1,aux,rhat,xlims;
    
    rhat=3;
    
    //Setting global 'myrhat_g' below
    //This global 'myrhat' will be referenced
    //in 'integrandMG'
    myrhat_g=rhat;
    
    xlims = { 1, 0 };
    z1 = intsimp(&integrandMG,xlims,1e-8);
    aux=z1^2;
    
    retp(aux);
endp;

proc integrandMG(t);
    local hy,aAlfa,bBeta,rhat;
    
    {aAlfa,bBeta}=alfabetaMG(pars,t);
    
    hy=exp(aAlfa+bBeta*myrhat_g);
    retp(hy);
endp;

proc (2)=alfabetaMG(pars,t);
    local aAlfa,bBeta,a,b,c,d,u;
    a=pars[1];
    b=pars[2];
    c=pars[3];
    d=pars[4];
    u=pars[5];
    aAlfa=a+b*t;
    bBeta=c+d*t;
    retp(aAlfa,bBeta);
endp;

proc startval();
    local a,b,c,d,u,pars;
    a=0.2;
    b=0.9;
    c=0.3;
    d=1;
    u=0;
    pars=a|b|c|d|u;
    retp(pars);
endp;

link

aptech

1,773

0

Yes that is perfect, thanks a lot.

However, I am afraid I have a final question.

In my final code, this example is part of an optimization procedure. So the complication that this brings is that the "pars" which now are global, will not be anymore and have to be passed on to each procedure as inputs.

The problem here is in the procedures that are the input of the integrals. There I cannot add any other input but "t" (my integrating variable) but inside these I need to call another procedure that gives me the function to be integrated. This function needs something more than "t" which it cannot get apparently and I get an error.

To put it in terms of the previous example, in

proc integrandMG(t);
    local hy,aAlfa,bBeta,rhat;
    
    {aAlfa,bBeta}=alfabetaMG(pars,t);
    
    hy=exp(aAlfa+bBeta*myrhat_g);
    retp(hy);
endp;

I get the error of "pars" is not defined. This is because it is not global anymore.

Any suggestions on how to solve this?

Thanks.

link

Babi

6