function [Call Put]= CEVprice(So,strike,r,T,sigma,t,Y,q)
v=sigma^2/(2*(r-q)*(Y-1)).*(exp(2*(r-q)*(Y-1).*(T-t))-1);
b=1/(1-Y);
warning('off','MATLAB:divideByZero')
% The CEV vanilla option pricing formulae in (Hull 2003),
% don't work when Y=0. I.e. we end up dividing by zero in
% the formulae. The "if" statement below, is taking into
% account the limiting values of the option pricing formulae
% in this case. The same problem occurs when T=0.
% The limiting value of the option pricing formulae in this
% case are also taken into account in the "if" statements
% below. Hence, though a Warning occurs, one need not worry
% about it - for this reason it has been switched off.
if (0<Y)&(Y<1)
a=(strike.*exp(-(r-q).*(T-t))).^(2*(1-Y))./((1-Y).^2*v);
c=So^(2*(1-Y))./((1-Y)^2*v);
a(find(isinf(a)))=0;
c(find(isinf(c)))=inf;
Call=So.*exp(-q.*(T-t)).*(1-ncx2cdf(a,b+2,c))-...
strike.*exp(-r.*(T-t)).*ncx2cdf(c,b,a);
Call=(Call>0).*Call;
a=(strike.*exp(-(r-q).*(T-t))).^(2*(1-Y))./((1-Y).^2*v);
c=So^(2*(1-Y))./((1-Y)^2*v);
a(find(isinf(a)))=inf;
c(find(isinf(c)))=0;
Put=strike.*exp(-r.*(T-t)).*(1-ncx2cdf(c,b,a))-...
So.*exp(-q.*(T-t)).*ncx2cdf(a,b+2,c);
Put=(Put>0).*Put;
elseif (Y>1)
a=(strike.*exp(-(r-q).*(T-t))).^(2*(1-Y))./((1-Y).^2*v);
c=So^(2*(1-Y))./((1-Y)^2*v);
a(find(isinf(a)))=inf;
c(find(isinf(c)))=0;
Call=So.*exp(-q.*(T-t)).*(1-ncx2cdf(c,-b,a))-...
strike.*exp(-r.*(T-t)).*ncx2cdf(a,2-b,c);
Call=(Call>0).*Call;
a=(strike.*exp(-(r-q).*(T-t))).^(2*(1-Y))./((1-Y).^2*v);
c=So^(2*(1-Y))./((1-Y)^2*v);
a(find(isinf(a)))=0;
c(find(isinf(c)))=inf;
Put=strike.*exp(-r.*(T-t)).*(1-ncx2cdf(a,2-b,c))-...
So.*exp(-q.*(T-t)).*ncx2cdf(c,-b,a);
Put=(Put>0).*Put;
elseif Y==1
[Call,Put] = blsprice(So, strike, r, T-t, sigma,q);
else
'error'
end