[x,coef,sol]=leuclag1(n,bas,Ps,u0,alpha,beta,gamma)

function [x,coef,sol]=leuclag1(n,bas,Ps,u0,alpha,beta,gamma)

[x,w,T,iT,X,D]=lagtool(n,bas);

cinit=zeros(3*n,1);cinit(1)=alpha;cinit(n+1)=beta;cinit(2*n+1)=gamma;
options=optimset('Display','iter','TolX',1.e-12,'TolFun',1.e-12,'MaxFunEvals',100000);
[coef,fval] = fsolve(@f,cinit,options);

sol(:,1)=T*coef(1:n)+u0(1);sol(:,2)=T*coef(n+1:2*n)+u0(2);sol(:,3)=T*coef(2*n+1:3*n)+u0(3);

function F=f(coef)
c1=T*coef(1:n)+u0(1);c2=T*coef(n+1:2*n)+u0(2);c3=T*coef(2*n+1:3*n)+u0(3);
P=param;
a=P(1,1);b=P(1,2);c=P(1,3);
A=P(2,1);B=P(2,2);C=P(2,3);
g=P(3,1);G=P(3,2);h=P(3,3);

f1=(a*(c1+c2)./(1+b*(c1+c2+c3))./(c1+c2+g*c3)-c).*c1;
f2=(A*(c1+c2)./(1+B*(c1+c2+c3))./(c1+c2+G*c3)-C).*c2;
f3=(a*c3./(1+b*(c1+c2+c3))./(c3+h*(c1+c2))-c).*c3;

F(1:n)=D*coef(1:n)-iT*f1;%F(n)=sum(coef(1:n))-alpha;
F(n+1:2*n)=D*coef(n+1:2*n)-iT*f2;%F(2*n)=sum(coef(n+1:2*n))-beta;
F(2*n+1:3*n)=D*coef(2*n+1:3*n)-iT*f3;%F(3*n)=sum(coef(2*n+1:3*n))-gamma;
z=Ps*[sum(coef(1:n));sum(coef(n+1:2*n));sum(coef(2*n+1:3*n))]-[alpha;beta;gamma];
F(n)=z(1);F(2*n)=z(2);F(3*n)=z(3);

end
end