Solving DAE

I'm having a trouble to solve the following system of DAE syetem when n<0.5 here 0<n<=1 and P is a given parameter.

dw/ds+(n/(n+1)) k^((1-n)) k^3-(alpha*k+p)=0

dv/ds=w

0 dk/ds=v k^((1-n))-k (equivalently v = k|k|^(1-n)-k)

dtheta/ds = k

dx/ds = cos(theta)

dy/ds= sin(theta) where alpha is a constant to be determined, with the BCs theta(0)=0, theta(2*pi/N) = 2*pi/N , k'(0)=0, k'(2*pi/N) = 0, x(0)=y(0)=0.

I solve this system by converting it to an initial value problem as below theta(0)=0; k(0)=beta; v(0) = beta^n if beta>=0 and -(-beta)^n if beta<0 w(0)=0, x(0)=y(0)=0;

I solve this initial value problem (as a DAE because of the second equation) by odesolver15s (to handle DAE) for given guesses of alpha and beta, while n is fixed and P needs to be varying then I use Newtons method to update the values of alpha and beta to satisfy the BCs theta(2*pi/N) = 2*pi/N and w(2*pi/N)=0. As far as I know BVP solver in Matlab doesn’t have the capability to solve a DAE system. I get initial guesses of alpha and beta with some other method and increase the value of P and get new alpha beta and keep moving up to certain pressure. My code is working fine for some initial P values but the solver starts complaining (Warning: Failure at t=3.053276e+000. Unable to meet integration tolerances without reducing the step size below the smallest value allowed (7.105427e-015) at time t. ). The problem is due to the situation that k(s) become zero if you look at the figures produce by my code. I need someone’s help if he/she able to rectify this problem in the code or in mathematical equations what I'm using. I'm thankful in advance for someone who is gonna help me in this regard. My code is here as well to look at it

function DAEquation

global k N k=.3;% k could be 0<k<=1 N=2; %P=3; format long double('precision') delta=.001; %stepsize for Newton's method col=1; %initial guesses alpha= -0.647821786142291; beta = 1.698977950392436; y2not=[]; kmin = []; pre = [.93,1.27,1.8,3];

m1=length(pre); for i1=1:m1 press= pre(i1):.01:pre(i1+1); n1=length(press);

for n2=1:n1 P= press(n2); % Start of Newton's method %initialize parameters for Newton's method alpha1=0; beta1=0; var = [alpha; beta]; while (abs(var(1) -alpha1)+ abs(var(2)-beta1))>1e-3

alpha =var(1); beta=var(2); % to approximate Jacobian by second order approximation for p=1:2 if p==1 alpha1=alpha ; beta1=beta; alpham=alpha+delta; betam = beta; end if p==2 alpha1=alpha ; beta1=beta; alpham=alpha-delta; betam = beta;

end for q=1:2 y(1)=0; y(2)=betam;

if betam>=0 y(3)=betam^k; elseif betam<0 y(3)=-(-betam)^k; end y(4)=0; y(5)=0; y(6)=0;

tspan=[0,2*pi/N]; % Partition of time domain %---------------------------------------- % Generating mass matrix %---------------------------------------- M = eye(6); M(2,2)= 0; options = odeset('Mass',M,'MStateDependence','none','MassSingular','yes','reltol',1e-7,'abstol',1e-7); [s,theta] = ode15s(@system3var,tspan,y,options,alpham,betam,P); %Calling the odesolver 15s

y1 = theta(:,1); a(p,q)=y1(end);

          y4 = theta(:,4);
          b(p,q)=y4(end);

      if p==1
          betam=beta+delta;
          alpham = alpha1;
      elseif p==2
          betam=beta-delta;
          alpham = alpha1;
      end

      end
      %alpha=alpha+delta;
  end
  %Jacobian matrix elements

J(1,1) = -(a(2,1)-a(1,1))/delta/2; J(1,2) = (a(1,2)-a(2,2))/delta/2; J(2,1) = -(b(2,1)-b(1,1))/delta/2; J(2,2) = (b(1,2)-b(2,2))/delta/2;

var = [alpha1; beta1];

y(1)=0; y(2)=beta1; %y(3)=beta*(abs(beta))^(k-1); if beta1>=0 y(3)=beta1^k; elseif beta1<0 y(3)=-(-beta1)^k; end y(4)=0; y(5)=0; y(6)=0;

tspan=[0,2*pi/N]; % Partition of time domain %finding solution at alpha and beta for Newtons method %---------------------------------------- % Generating mass matrix %---------------------------------------- M = eye(6); M(2,2)= 0; options = odeset('Mass',M,'MStateDependence','none','MassSingular','yes','reltol',1e-5,'abstol',1e-5,'refine',16); [s1,theta1] = ode15s(@system3var,tspan,y,options,alpha1,beta1,P); %Calling the odesolver 15s

y11 = theta1(:,1); c=y11(end);

          y41 = theta1(:,4);
          d=y41(end);

fun = [c-2*pi/N; d]; jac=inv(J); var = var- inv(J)*fun; end P alpha=var(1) beta=var(2)

end %Plotting theta and s figure(1) plot(s,theta(:,1)); hold on; %ploting k vs s figure(2) plot(s,theta(:,2)); hold on;

figure(4) plot(s,theta(:,3)); hold on; figure(5) plot(s,theta(:,4)); hold on;

figure(3) if col==1 plot([-theta(:,5),theta(:,5)],[theta(:,6),theta(:,6)],'m'); hold on; axis equal; end if col==2 plot([-theta(:,5),theta(:,5)],[theta(:,6),theta(:,6)],'r'); hold on; axis equal; end if col==3 plot([-theta(:,5),theta(:,5)],[theta(:,6),theta(:,6)],'b'); hold on; axis equal; end

if col==4 plot([-theta(:,5),theta(:,5)],[theta(:,6),theta(:,6)],'k'); hold on; axis equal; end

if col>4 plot([-theta(:,5),theta(:,5)],[theta(:,6),theta(:,6)],'c'); hold on; axis equal; end if col>10 plot([-theta(:,5),theta(:,5)],[theta(:,6),theta(:,6)],'y'); hold on; axis equal; end

col=col+1;

end end

function ode = system3var(t,y,alpha,beta,P) global k; if y(2)>=0 z=(y(2))^(k); elseif y(2)<0 z=-(-y(2))^(k); end ode(1) = y(2); ode(2) = y(3)- z; ode(3) = y(4); ode(4) = (P + alpha*y(2) - (k/(k+1))*z*(abs(y(2)))^(2)); ode(5) = cos(y(1)); ode(6) = sin(y(1));

ode=ode'; end