dy=evolution(time,y,strainrate,a,b,A0,H0,K,G,GnH,T)

function dy=evolution(time,y,strainrate,a,b,A0,H0,K,G,GnH,T)
%This is the evolution equation for Be, B, J, and L. The pseudo-vector y is composed of Be33, Be11, J, and L11, 
%and the pseudo-vector out contains the time derivatives of Be33, Be11, and J, and the algebraic dynamic boundary condition T11==0.
%Units: tau0[MPa], Vact[nm3]

%This is the true strain vector
truestrain=strainrate * time;

%In this model, the activation volume is a function of invariants of the strain tensor. 
%To edit the exact functional form of the activation volume, please edit the file Vact.m
Vactivation=Vact(a,b,truestrain);

tau0=1.380658 * 1e-2 * T / Vactivation; 

Be=[y(2) 0 0; 0 y(2) 0; 0 0 y(1)];
J=y(3);
L=[y(4) 0 0; 0 y(4) 0; 0 0 strainrate];

%Evolution equation for Be and J:
dBedt=(dev(symm(L))-dplas(sigma(J,Be,K,G),tau0,A0,H0,T))*Be+Be*(dev(symm(L))-dplas(sigma(J,Be,K,G),tau0,A0,H0,T));
dJdt=J*trace(symm(L));

%calculating the deviatoric part of B
devB=dev(B(truestrain,0.5));

%Boundary condition in uniaxial deformation: 0==T(1,1). 
%Note: Here we calculate the right side of the boundary condition, which is set to 0 in WTS.m via the matrix M.
T11=K*(J-1)+ (1/3) * G * (Be(1,1)-Be(3,3)) + GnH * devB(1,1);

%Pseudo-vector y out. dy(4) contains the right side of the boundary condition 0==T(1,1).
dy=[dBedt(3,3);dBedt(1,1);dJdt;T11];