%{
Function equations returns the vector value of the Euler Lagrange Equations
The function takes as parameter a solution array u, the point of interest i,j, 
with grid spacing dx, and the material parameters a,kappa2,kappa3 and 
interpolants Lam1, Lam2, and Lam3.
%}

function v = equationsMS(u,i,j,dx,a,kappa2,kappa3,Lam1,Lam2,Lam3)
    v = zeros(3,1); 
    %Full RHS of equation
    v(1) = -Lam1(u(1,i,j),u(2,i,j),u(3,i,j)) + (4/3)*a*u(1,i,j) -(2/9)*kappa3*(yderiv(u,1,i,j,dx)^2 + 3*yderiv(u,2,i,j,dx) + 3*yderiv(u,3,i,j,dx)^2 - 2*(xderiv(u,1,i,j,dx)^2 + 3*(xderiv(u,2,i,j,dx)^2 + xderiv(u,3,i,j,dx)^2))) + (4/3)*(xxderiv(u,1,i,j,dx) + yyderiv(u,1,i,j,dx)) + (2/3)*kappa2*((4/3)*xxderiv(u,1,i,j,dx) + (1/3)*yyderiv(u,1,i,j,dx) - yyderiv(u,2,i,j,dx) + xyderiv(u,3,i,j,dx)) + (2/9)*kappa3*(4*xderiv(u,1,i,j,dx)^2 - 2*yderiv(u,1,i,j,dx)^2 + 6*yderiv(u,2,i,j,dx)*yderiv(u,1,i,j,dx) + 6*yderiv(u,3,i,j,dx)*xderiv(u,1,i,j,dx) + u(1,i,j)*(4*xxderiv(u,1,i,j,dx) - 2*yyderiv(u,1,i,j,dx)) + 6*u(2,i,j)*yyderiv(u,1,i,j,dx) + 12*u(3,i,j)*xyderiv(u,1,i,j,dx));
    v(2) = -Lam2(u(1,i,j),u(2,i,j),u(3,i,j)) + 4*a*u(2,i,j) - (2/3)*kappa3*(yderiv(u,1,i,j,dx)^2 + 3*(yderiv(u,2,i,j,dx)^2 + yderiv(u,3,i,j,dx)^2))+ 4*(xxderiv(u,2,i,j,dx) + yyderiv(u,2,i,j,dx)) - 2*kappa2*((1/3)*yyderiv(u,1,i,j,dx) - (yyderiv(u,2,i,j,dx) + xyderiv(u,3,i,j,dx))) + (2/9)*kappa3*(12*xderiv(u,1,i,j,dx)*xderiv(u,2,i,j,dx) - 6*yderiv(u,1,i,j,dx)*yderiv(u,2,i,j,dx) + 18*yderiv(u,2,i,j,dx)^2 + 18*xderiv(u,3,i,j,dx)*yderiv(u,2,i,j,dx) + 6*u(1,i,j)*(2*xxderiv(u,2,i,j,dx) - yyderiv(u,2,i,j,dx)) + 18*u(2,i,j)*yyderiv(u,2,i,j,dx) + 36*u(3,i,j)*xyderiv(u,2,i,j,dx));
    v(3) = -Lam3(u(1,i,j),u(2,i,j),u(3,i,j)) + 4*a*u(3,i,j) - (4/3)*kappa3*(yderiv(u,1,i,j,dx)*xderiv(u,1,i,j,dx) + 3*yderiv(u,2,i,j,dx)*xderiv(u,2,i,j,dx) + 3*yderiv(u,3,i,j,dx)*xderiv(u,3,i,j,dx)) + (4 + 2*kappa2)*(xxderiv(u,3,i,j,dx) + yyderiv(u,3,i,j,dx)) + kappa2*((2/3)*xyderiv(u,1,i,j,dx) + 2*xyderiv(u,2,i,j,dx)) + (2/9)*kappa3*(12*xderiv(u,1,i,j,dx)*xderiv(u,3,i,j,dx) - 6*yderiv(u,1,i,j,dx)*yderiv(u,3,i,j,dx) + 18*yderiv(u,2,i,j,dx)*yderiv(u,3,i,j,dx) + 36*xderiv(u,3,i,j,dx)*yderiv(u,3,i,j,dx) + 6*u(1,i,j)*(2*xxderiv(u,3,i,j,dx) - yyderiv(u,3,i,j,dx)) + 18*u(2,i,j)*yyderiv(u,3,i,j,dx) + 36*u(3,i,j)*xyderiv(u,3,i,j,dx));
end