function [alphan,An,Bn] = Mandel_Newton 
% Zhuoran Wang, SYSU; James Liu, ColoState; 2017/07--2021/04 

%% Number of roots
terms = 1000;

%% Parameters for Mandel's problem 
% E = 9*K/(3*K+mu)*mu;  % Young's modulus
% nu = (3*K-2*mu)/(2*(3*K+mu));  % Poisson's coefficient
nu = 0.4;  E = 1e4;  
lambda = E*nu/(1-2*nu)/(1+nu);  mu = E/2/(1+nu); 

%% For Darcy 
kappa = 1e-6;  % Choices: 1e-6, 1e-3, 1, 1e3 
c0 = 0; 

%% More parameters 
% K = lambda + 2/3*mu; % Skeleton bulk modulus
% K = 1e-06;
% Ku = K+alpha^2/c0; % undrained version
% cf = kappa/c0*(K+4/3*mu)/(Ku+4/3*mu); % Fluid diffusivity coeff.
cf = kappa*(lambda + 2*mu);
% B = alpha/(c0*Ku); % Skempton's coeff.
% B = 3/(alpha*(1-2*nu))*(nuu-nu)/(1+2*nuu);
B = 1; 
% nuu = ((B-2*B*nu)/3+nu)/(1-2*B*(1-2*nu)/3);
nuu = (3*nu+B*(1-2*nu))/(3-B*(1-2*nu));

%% For coupling 
alpha = 1;  

%% More parameters 
a = 1;  b = 1;  % Domain size 
F = 1;  
T = 50;  % Final time 

%% Newton method
 f = @(x) tan(x)-(1-nu)/(nuu-nu)*x;
fp = @(x) sec(x)^2-(1-nu)/(nuu-nu);
for m = 1:terms
  x0 = (2*m-1)*pi/2 - 1E-5;  % Initial guess
  N = 1000; 
  tol = 1E-16;
  x(1) = x0;  % Set initial guess
  n = 2; 
  nfinal = N + 1; 
  while (n <= N + 1)
    fe = f(x(n - 1));
    fpe = fp(x(n - 1));
    x(n) = x(n - 1) - fe/fpe;
    if (abs(fe) <= tol)
      nfinal = n; 
      break;
    end
    n = n + 1;
  end
  alphan(m) = x(nfinal);
end

An = sin(alphan)./(alphan - sin(alphan).*cos(alphan));
Bn = cos(alphan)./(alphan - sin(alphan).*cos(alphan));
    
% % Graph exact pressure p and sigma_yy    
x = linspace(0,1); % change to h segmants.
NumPts = numel(x); % also change segmants.
% t = linspace(0,100,1000);
t = [0.1,1,10,50];
% t = [0.01,0.1,1,5]
size_t = size(t,2);
for p = 1:size_t
    for q = 1:NumPts
        for i=1:NumPts  % modify
             xhat = x(i)/a; that = cf*t(p)/(a^2);
             for j = 1:terms
                     px(i,j) = An(j)*exp(-alphan(j)^2*that)* ...
                             (cos(alphan(j)*xhat) - cos(alphan(j)));
                     sigma(i,j) = An(j)*exp(-alphan(j)^2*that)* ...
                     (cos(alphan(j)) - B*(nuu-nu)/(1-nu)*cos(alphan(j)*xhat));    
             end
        end
        phat_sub = 2*sum(px,2);
        sigmahat_yy_sub = -1+2*sum(sigma,2);
    end
    phat(p,:) = phat_sub;
    sigmahat_yy(p,:) = F/a*sigmahat_yy_sub;
end
% phat1 = 2*sum(px,2);
% sigmahat_yy = -1+2*sum(sigma,2);
% plot(x,phat1,'red');
% hold on;
% 
plot(x,phat(1,:),'black');
hold on;
plot(x,phat(2,:),'red');
hold on;
plot(x,phat(3,:),'blue');
hold on;
plot(x,phat(4,:),'green');
hold on;
% 
% figure(2);
% plot(x,sigmahat_yy(4,:),'black');
% hold on;

% % scatter(x,phat1,'*');
% hold on;
% % size(phat)
% s = mesh(x,t,phat) % normalized pressure % mesh
% s.FaceColor = 'flat';

end