function [Dcons,sig,epsE,epsP,sigmay,q,itnum] =Von_Mises_Isotropic_hardening(G,lamda,sigmay,h,epsEtr,epsP)
maxit=25; tol=1e-15;itnum=0;
eye1=[1 1 1 0 0 0]';eyeS=[eye(3,3) zeros(3,3);zeros(3,3) eye(3,3)*2];eyeD=eyeS-eye1*eye1'/3;
sigmay0=sigmay;epsE=epsEtr;De=eye(6,6)*2*G+eye1*eye1'*lamda;Dcons=De;
sig=De*epsEtr;s=(eye(6,6)-eye1*eye1'/3)*sig;
s_mat=[s(1) s(4) s(6);s(4) s(2) s(5);s(6) s(5) s(3)];J2=simplify(trace(s_mat*s_mat)/2);q=sqrt(3*J2);
f=q-sigmay; 
if (f>tol)
        res=zeros(7,1); res(8)=f; itnum=0; dlamda=0;dalpha=0;
        dJ2sig=eyeS*s;ddJ2sig=eyeD;eN=3/2/q*dJ2sig;deN=3/2*ddJ2sig/q-9/4/q^3*dJ2sig*dJ2sig';
        while (itnum<maxit) && ((norm(res(1:6))>tol) || (abs(res(7))>tol) || (abs(res(8))>tol))
            A=[eye(6,6)+dlamda*deN*De zeros(6,1) eN;zeros(1,6) 1 -1;eN.'*De -h 0];
            dx=-inv(A)*res; epsE=epsE+dx(1:6); dalpha=dalpha+dx(7);dlamda=dlamda+dx(8); 
            sigmay=h*dalpha+sigmay0;sig=De*epsE;s=(eye(6,6)-eye1*eye1'/3)*sig;
            s_mat=[s(1) s(4) s(6);s(4) s(2) s(5);s(6) s(5) s(3)];J2=simplify(trace(s_mat*s_mat)/2);q=sqrt(3*J2);
            dJ2sig=eyeS*s;ddJ2sig=eyeD;eN=3/2/q*dJ2sig;deN=3/2*ddJ2sig/q-9/4/q^3*dJ2sig*dJ2sig';
            res=[epsE-epsEtr+dlamda*eN;dalpha-dlamda; q-sigmay];
            itnum=itnum+1;
        end
        depsP=dlamda*eN;epsP=epsP+depsP;
        B=inv(A);Dcons=B(1:6,1:6);
end