%{
This function solves the 2D nematic tactoid problem
for grid size M, topological charge of boundary conditions m, 
initial radius of tactoid or defect Rin,
system radius of Rout, anisotropy constants kappa2 and kappa3, 
temperature parameter a, SOR parameter dt,
convergence factor err, maximum number of iterations numits, and 
MATLAB interpolants Lam1, Lam2, Lam3.
This version of the function includes up to cubic order in the
gradient free energy to break the Frank degeneracy.
%}
function [u,E] = tactoidSolnCubeMS(M,m,Rin,Rout,kappa2,kappa3,a,Sn,dt,err,numits,Lam1,Lam2,Lam3)
    dx = (2*Rout)/M; %grid spacing
    x = -Rout:dx:Rout; %x vector for system
    y = x; %use a square

    %Initialize u
    S = zeros(M+1,M+1);
    S(:,:) = Sn;
    
    %Check if bulk is nematic or coexistence
    if a >= 3.41
        for i=2:M
            for j=2:M
		if sqrt(x(i)^2 + y(j)^2) < Rout - 0.001
	            S(i,j) = Sn*(1 - exp(-sqrt(x(i)^2 + y(j)^2)/Rin));
                end
            end
	end
    else
        for i=2:M
	    for j=2:M
	        S(i,j) = (Sn/2)*(1 - tanh((Rin - sqrt(x(i)^2 + x(j)^2))/10));
	    end
        end
    end
    
    
    P = zeros(M+1,M+1);

    phi = zeros(M+1,M+1);
    
    for i=1:M+1
	for j=1:M+1
	    phi(i,j) = m*atan2(y(j),x(i));
	end
    end
    
    u = zeros(3,M+1,M+1);
    for i=1:M+1
	for j=1:M+1
	    u(1,i,j) = S(i,j) - (3/2)*(S(i,j) - P(i,j))*sin(phi(i,j))^2;
	    u(2,i,j) = P(i,j) + (1/2)*(S(i,j) - P(i,j))*sin(phi(i,j))^2;
	    u(3,i,j) = (1/2)*(S(i,j) - P(i,j))*sin(2*phi(i,j));
	end
    end


    %Call function getEnergy2 to get the energy
    E0 = getTacEnergyMS(x,u,a,kappa2,kappa3);
    E = E0;
    
    %Big loop, runs numits times
    for inc1 = 1:numits
        for i=2:M
            for j=2:M
		if sqrt(x(i)^2 + y(j)^2) < Rout - 0.001
	    	    u(:,i,j) = u(:,i,j) - dt*(equationsMS(u,i,j,dx,a,kappa2,kappa3,Lam1,Lam2,Lam3) ./ jacobian2MS(u,i,j,dx,a,kappa2,kappa3));
		end 
            end
	end
        %Calculate Energy
        E = getTacEnergyMS(x,u,a,kappa2,kappa3);
        if abs(E0 - E) < err*(abs(E0) + abs(E) + eps) %If energy stops changing, stop iterating
            disp(['Convergence after ',num2str(inc1),' iterations']);
            break;
        end
	E0 = E;

    end
    if numits == 0
	E = E0;
    end
    
end