%solve ODE 
[t,x] = ode23s('p426sbfun',[0;10000],[1544.0,0.0,0.0,0.0]); 
%plot(t,x(:,1),'o',t,x(:,2),'d',t,x(:,3),'t'); 
na=x(:,1);
nb=x(:,2);
nc=x(:,3);
nd=x(:,4);
Vi=200;
tind=length(t);
for i=1:tind,
    V=Vi+0.05*t(i);
end
V=V';
ca=na./V;
cb=nb./V;
cc=nc./V;
cd=nd./V;
xa=(1544-na)/1544;