% read input parameters
lleq_in

if ~restart
    
clear xy
clear ds
%generate grid and del_star matrix
SIgrid
delstar

%LU decomposition of del_star matrix
[dsL,dsU]=lu(ds);

if ~remeshed
    opts.disp=0;
    
    % calculate Taylor state (eigenvector of del_star)
    [psi,mlambda2]=eigs(ds,1,-121,opts);
    
    lambda=sqrt(-mlambda2);
    lambda_bar=lambda;
end

end

% if Taylor state is desired, this is the solution
psi=abs(psi)/(max(abs(psi)));

%otherwise iterate solution with Taylor state as starting guess
if ~(alpha==0)
        
for iter=1:max_iter 

% left hand side

old_psi=psi;

lhs=-((1+alpha.*(psi-1))+beta.*(psi.^2-1)).*(1+alpha.*(2.*psi-1)+beta.*(3*psi.^2-1)).*psi;


y=dsL\lhs;
psi=dsU\y;

lambda_bar=1/sqrt(max(abs(psi)));
psi=abs(psi)/(max(abs(psi)));


error=sqrt(sum((psi-old_psi).^2)/(length(psi)));
rms_err(iter)=error;

if rms_err(iter)<rms_tol
break
end

end
psi=abs(psi)/(max(abs(psi)));
iterations=iter
rms_error=error

end

%put solution onto the grid and plot
eq_block
flux_plot