clear all
close all

m=[0.3,0.5,0.7];
subf=m.^(1/2);
alpha=-22.2
delta=10^4;
tao=1.67*10^(-4);
rho=10;
A=(-2)*alpha*tao;

N=100;

funstr='1/sqrt((rho+w-delta/(rho+1))*(1-y^2)+rho/2*(y^4-1)-w*(1/y^2-1)-delta/rho*ln((rho*y^2+1)/(rho+1)))';

for ii=1:length(m)
w=rho*m(ii)/2-delta*m(ii)/(rho+1)/(1-m(ii))-delta*m(ii)/rho/(m(ii)-1)^2*log((1+rho*m(ii))/(1+rho));

y1=linspace(subf(ii),1,N);

fun=subs(funstr,{'w','rho','delta'},{w,rho,delta})

for jj=1:N

s1(jj)=quadl(inline(vectorize(char(fun))),subf(ii),y1(jj));
s1(jj)=1/sqrt(A)*s1(jj);
end
ss(ii,:)=[fliplr(-s1) s1];
yy(ii,:)=[fliplr(y1.^2) y1.^2];

end

disp('------------计算结束---------------')
%%画图
figure
hold on

axis([-10 10 0 1.1]);
line(ss(1,:),yy(1,:)*30/4,'color','k');
line(ss(2,:),yy(2,:)*30/4,'color','k','linestyle',':');
line(ss(3,:),yy(3,:)*30/4,'color','k','linestyle','-.');
line(ss(1,:),yy(1,:)*10/4,'color','k','color','k');
line(ss(2,:),yy(2,:)*10/4,'color','k','linestyle',':');
line(ss(3,:),yy(3,:)*10/4,'color','k','linestyle','-.');

legend('m=0.3','m=0.5','m=0.7');