[T33,truestrain,J,S33,R33]=WTS(K,G,GnH,a,b,A0,H0,straindata,strainrate,T,tol)

function [T33,truestrain,J,S33,R33]=WTS(K,G,GnH,a,b,A0,H0,straindata,strainrate,T,tol)
%This function determines the total true axial stress in a glassy polymer in uniaxial deformation according to the WTS model

% A constant, singular mass matrix for uniaxial deformation satisfying the DAE M*y'=f(y,t). 
M = [1 0 0 0; 0 1 0 0; 0 0 1 0; 0 0 0 0];

% Initial conditions at t=0 of the pseudo vector y=[Be(3,3); Be(1,1); J; T(1,1)]. 
y0=[1; 1; 1; 0];

%time intervall
tspan=straindata/strainrate;

%solver options
options=odeset('Mass',M,'AbsTol',[1e-5 1e-5 1e-5 tol],'MassSingular','yes','MStateDependence','none','InitialSlope',evolution(0,y0,strainrate,a,b,A0,H0,K,G,GnH,T));

%start solving the pseudovector Y from the DAE for the viscoelastic stress response
[Time,Y]=ode15s(@evolution,tspan,y0,options,strainrate,a,b,A0,H0,K,G,GnH,T);

%obtaining the axial stress, strain, and relative volume change J for uniaxial deformation
S33=K*(Y(:,3)-1)+ (2/3)*G*(Y(:,1)-Y(:,2));
truestrain=strainrate * Time;
J=Y(:,3);

%adding a neo-Hookean stress response with modulus GnH to account for the elastic contribution to strain hardening
%Note: GnH has no temperature dependence in this code
for tel=1:length(truestrain)
    devB=dev(B(truestrain(tel),0.5));
    R33(tel,1)=GnH * devB(3,3);
end

%total stress T=S+R
T33=S33 + R33;