%{
Function runs a backward Euler method to solve the one dimensional 
solution for the Maier-Saupe model. Inputs are number of discrete
spaces M, system size Rout, temperature paramter a,
nematic order parameter Sn, angle of director phi, anisotropy constant
kappa2, convergence criteria err, max iterations err, and 
MATLAB interpolants Lam1, Lam2, Lam3
%}
function [eta,mu,nu,E] = backEulerMS1D(M,Rout,a,Sn,phi,kappa2,err,numits,Lam1,Lam2,Lam3)
    %Initialize values
    x = linspace(-Rout,Rout,M);
    dx = x(2)-x(1);
    dt = 0.1*dx^2;
    q = dt/(dx^2);
    
    %Build initial condition
    S = 0.5*Sn*(1 - tanh(x/10));
    S(1) = Sn;
    S(M) = 0;
    
    P = zeros(1,M);
    
    eta0 = S - 1.5*(S - P)*sin(phi)^2;
    mu0 = P + 0.5*(S - P)*sin(phi)^2;
    nu0 = 0.5*(S - P)*sin(2*phi);
    
    %Dirichlet Boundary conditions
    A = eta0(1);
    B = mu0(1);
    C = nu0(1);
    
    E0 = getE1D(x,eta0,mu0,nu0,a,kappa2);
    eta = zeros(1,M);
    mu = zeros(1,M);
    nu = zeros(1,M);
    for its = 1:numits
        %Build tridiagonal systems
        aeta = zeros(1,M);
        beta = zeros(1,M);
        ceta = zeros(1,M);
        reta = zeros(1,M);
        beta(1) = 1;
        beta(M) = 1;
        reta(1) = A;
        amu = zeros(1,M);
        bmu = zeros(1,M);
        cmu = zeros(1,M);
        rmu = zeros(1,M);
        bmu(1) = 1;
        bmu(M) = 1;
        rmu(1) = B;
        anu = zeros(1,M);
        bnu = zeros(1,M);
        cnu = zeros(1,M);
        rnu = zeros(1,M);
        bnu(1) = 1;
        bnu(M) = 1;
        rnu(1) = C;
        for j=2:M-1
            aeta(j) = -q*((4/3) + (8/9)*kappa2);
            beta(j) = 1 + q*((8/3) + (16/9)*kappa2);
            ceta(j) = -q*((4/3) + (8/9)*kappa2);
            reta(j) = eta0(j) + dt*(-Lam1(eta0(j),mu0(j),nu0(j)) + (4/3)*a*eta0(j));
            amu(j) = -q*4;
            bmu(j) = 1 + q*8;
            cmu(j) = -q*4;
            rmu(j) = mu0(j) + dt*(-Lam2(eta0(j),mu0(j),nu0(j)) + 4*a*mu0(j));
            anu(j) = -q*(4 + 2*kappa2);
            bnu(j) = 1 + q*(8 + 4*kappa2);
            cnu(j) = -q*(4 + 2*kappa2);
            rnu(j) = nu0(j) + dt*(-Lam3(eta0(j),mu0(j),nu0(j)) + 4*a*nu0(j));
        end
        %Solve tridiagonal system for next iteration
        eta = tridag(aeta,beta,ceta,reta);
        mu = tridag(amu,bmu,cmu,rmu);
        nu = tridag(anu,bnu,cnu,rnu);
        E = getE1D(x,eta,mu,nu,a,kappa2);
        if abs(E-E0) < err*(abs(E) + abs(E0) + eps)
            disp(['Convergence after ',num2str(its),' iterations']);
	    return;
        end
        E0 = E;
        eta0 = eta;
        mu0 = mu;
        nu0 = nu;
    end
end