function s=simcev(n,r,sigma,So,T,gam)
% simulates n sample paths of a CEV process on the interval [0,T] all with
% the same starting value So. assume gamma != 1.
%  s=simcev(10,.05,.2,10,1,.8);
Yt=ones(n,1)*(So^(1-gam))/(1-gam);  y=Yt;
dt=T/1000; c1=r*(1-gam)/sigma; c2=gam*sigma/(2*(1-gam));
dw=normrnd(0,sqrt(dt),n,1000);
for i=1:1000
    v=find(Yt);   % selects positive components of Yt for update
    Yt=max(0,Yt(v)+(c1.*Yt(v)-c2./Yt(v))*dt+dw(v,i));
    y=[y Yt];
end
s=((1-gam)*max(y,0)).^(1/(1-gam));  %transforms to St