function nup=nu_p()
   syms x_low x_up R sig s;    
   kappa=0.0858;
   p=0.5;
   h=hfun();
   g=gfun();
   h_x=hfun_deriv();
   g_x=gfun_deriv();
   nup=(kappa-1)*(1+(h_x/(1-p)*h*(g_x/g)-h_x));
end
function h=hfun()
   syms x_low x_up R sig s;    
   eta=0.27;
   mu=0.12;
   a=norm1(s);
   h=(1-mu-eta)*int(a,x_low,x_up)+R*int(a,x_up,inf)-x_low*int(a,eta*s,s)-mu*int(a,0,x_low);
end
function a=norm1(s)
a=normpdf(s);
end
function [h_x,h_sig]=hfun_deriv()
syms x_low x_up R sig s; 
eta=0.27;
ss=eta*s;
mu=0.12;
a=norm1(s);
b=norm2(x_up);
c=norm3(ss);
h_x=(1-mu-R-eta)*b;
h_sig=-x_low*(a-c);
end
function b=norm2(x_up)
b=normpdf(x_up);
end
function c=norm3(ss)
c=normpdf(ss);
end
function g=gfun()
   syms x_low x_up R sig s sp sr;    
   eta=0.27;
   p=0.5;
   sp=s.^(1-p);
   sr=(s-R).^(1-p);
   a=norm1(s);
   d=norm4(sp);
   e=norm5(sr);
   g=int(d,0,x_low)+x_low.^(1-p)*int(a,eta*s,s)+eta.^(1-p)*int(d,x_low,x_up)+int(e,x_up,inf);
end
function d=norm4(sp)
d=normpdf(sp);
end
function e=norm5(sr)
e=normpdf(sr);
end
function [g_x,g_sig]=gfun_deriv()
syms x_low x_up R sig s ss zz RR; 
eta=0.27;
p=0.5;
ss=eta*s;
zz=x_up.^(1-p);
RR=(x_up-R).^(1-p);
a=norm1(s);
f=norm6(zz);
j=norm7(ss);
i=norm8(RR);
g_x=eta.^(1-p)*f-i;
g_sig=x_low.^(1-p)*(a-j);
end
function f=norm6(zz)
f=normpdf(zz);
end
function j=norm7(ss)
j=normpdf(ss);
end
function i=norm8(RR)
i=normpdf(RR);
end