function Error=ShapePreserving(T,Y)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% ShapePreserving.m constructs the best convex-concave and otherwise free
%% constrained cubic splines intrepolating given data.
%
% Modified: 17 January, 2007
%
% Author: 
%
% Houduo Qi
% School of Mathematics
% University of Southampton
% Highfield
% Southampton SO17 1BJ
% England
%%
% Input
% T=[t_1, t_2,...,t_{N+2}]'; (column vector)
% Y=[y_1, y_2,...,y_{N+2}]'; (column vector)
%
% Output
% Error=norm of gradient of the objective function f
% figure: dashed line is the resulting shape-preserving cubic spline
%         solid line is the natural spline interpolating the given data

% References
%
% A. Dontchev, H.-D. Qi, L. Qi, and H. Yin:
% A Newton method for shape-prserving spline interpolation
% SIOPT, 13 (2): 588--602
%
% MATLAB code sptsl.m by John Burkardt was used to solve the tridiagonal
% linear equations appeared in Newton's method for the best convex case
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 [N,NC]=size(T);
 N=N-2;
  
%    Parameter settings
%
  rho=0.5; %(Armijo steplenth)
   sigma=0.1;
   delta=0;
   tol=1e-10;
%  kondmax=1e16; 
  tmin= 2^(-3);%1.e-16;

% Number of iteration 

   k=0;

% Number of line search used
  lsit=0;
% The maximum iteration allowed

  kmax=1000;

% 2nd divided differences

  D=zeros(N,1);

  for i=1:N
   Delta1=(Y(i+1)-Y(i))/(T(i+1)-T(i));
   Delta2=(Y(i+2)-Y(i+1))/(T(i+2)-T(i+1));   
   D(i)=1/(T(i+2)-T(i))*(Delta2-Delta1);
  end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% Data Removal Process Started (if you want to remove some data)
% % Data removal: Remove data Y(i+1) if D(i)==0
% % J1=zeros(N+2,1);
% 
% for i=1:(N+2)
%     J1(i)=i;
% end
% I=find(abs(D) <= 1.0e-5); % D(i) is thought to be zeros if |D(i)| <= small number
% 
% [s1, s2]=size(I);
% if s2==0 % I is empty
%    return;
% else
%    for i=1:s1
%        J1(I(i)+1)=0;
%    end
%    J=find(J1>0.5);
%    T=T(J);
%    Y=Y(J);
%    %%% Data removal process is completed
%    % Re-set other paratemetres associated with T
%    [N,NC]=size(T);
%    N=N-2;
%    % New 2nd divided differences
%   D=zeros(N,1);
% 
%   for i=1:N
%    Delta1=(Y(i+1)-Y(i))/(T(i+1)-T(i));
%    Delta2=(Y(i+2)-Y(i+1))/(T(i+2)-T(i+1));   
%    D(i)=1/(T(i+2)-T(i))*(Delta2-Delta1);
%   end
% 
% end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%End of data removal process%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

 
% Defining II, the classification vector

  II=zeros(N,1);
  
% Allow 0 values in D(i)
  if D(1)>0
    II(1)=1;
  elseif D(1)<0
    II(1)=-1;
  else
    II(1)=0;
  end

  for i=2:N
   if (D(i-1) > 0) & (D(i) >0)
     II(i)=1;
      elseif (D(i-1)<0) & (D(i) <0)
     II(i)=-1;
        else
     II(i)=0;
   end
  end
% STARTING POINT


  Lam=sign(D);   

% We are going to solve the equations:
% F(\lambda)=2d (see, [Andersson and Elfving [SIAM Scientifc 87, 1014]

% Calculate the gradient $F(Lam)$ and the function $f(\Lam)$ itself

   [F, f]=Ff(T, Lam, D, II);

  Fold=F;
  fold=f;
% Number of function evaluation

   nf=1;
   
% Norm of the gradient of f

   ngf=norm(Fold);

% Calculate the Jacobian matrix $M(\lambda)$

   [dM, offdM]=Jacob(T, Lam, II);

% Loop
%
  while (k < kmax) & (ngf > tol)
      k=k+1;
%%%%%%use sptsl to solve the tridiagonal linear equations
        dM=dM+min(delta, ngf); % diagonal vector of Mk
%     regularization with a small positive diagonal matrix
%     numerical performance shows that no regularization performs best
%     hence, in this code, delta=0.
        d=sptsl(N,dM,offdM,-Fold);
%%%%%%%% end of sptsl
%%%%%%line search%%%%%%%%%%%%%%%%%%%%%%%%%%%%
	t=1;
	Lamnew=Lam + t*d;
    [Fnew, fnew]=Ff(T, Lamnew, D, II);
	nf=nf+1;
    
    bound=sigma*Fold'*d;
    
       while (fnew > (fold + t*bound))  & (t>tmin) 
	      
	      lsit=lsit+1;
          t=rho*t;
	      Lamnew=Lam +t*d;
          [Fnew,fnew]=Ff(T, Lamnew, D, II);
	      
	      nf=nf+1;
	    end
%%%%%%%%% end of line search%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%Updates%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

	Lam=Lamnew;
	fold=fnew;
    Fold=Fnew;
	ngf=norm(Fold);
	[dM, offdM]=Jacob(T,Lam, II);
  end % end of while (first)

  Error=ngf;
% Output some information as the number of iterations
 fprintf('Number of iteration =%d\n', k);
 fprintf('Size of the problem =%d\n', N);
 fprintf('Number of function evaluations =%d\n', nf);
 fprintf('Number of line search called =%d\n', lsit);
 fprintf('Final value of the error =%10.1e\n', ngf);

% end
% Now generating the convex cubic spline
% The procedure is taken from [Andersson and Elfving]
% Let \bar \lambda=4\lambda/(t_{i+2)-t_i)
%
  Barlam=zeros(N,1);
  for i=1:N
      Barlam(i)=4*Lam(i)/(T(i+2)-T(i));
  end

 cubic=Cubic(T, Y, II, Barlam);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Enf of Main Program: ShapePreserving.m%%%%%%%%%%%%%%%%%%%%

 function [F,f]=Ff(T, Lam, D, II)

%The function fvalue uses one application of SIMPSON's rule to calculate
%numerically the integral which reaults in the functional value of $f(Lam)$
% BEND is the basic matrix on data T, that is the matrix BEND in the main
% program, BMID is a matrix on the data (T(i)+T(i+1))/2
% II is an indentification vector whose values are 1, -1, 0 and
% II(i)=1 means [t_i, t_{i+1}] belongs to $\Omega_1$
% II(i)=-1 means [t_i, t_{i+1}] belongs to $\Omega_2$
% II(i)=0 means [t_i, t_{i+1}] belongs to $\Omega_3$

% f(\lambda)=\int_a^b(\sum_{l=1}^N \lambda_l B_l)_+^2X_{\Omega_1} dt
%            +\int_a^b(\sum_{l=1}^N \lambda_l B_l)_-^2X_{\Omega_2} dt
%            +\int_a^b(\sum_{l=1}^N \lambda_l B_l)^2X_{\Omega_3} dt

 [N,NC]=size(T);
 N=N-2;
 f=0;
 
 % Use Gamma and Psi to calculate F (DQQY03: Eq.(6))
 
 % for II(i)==1 ([t_i, t_{i+1}] \in \Omgea_1
 % Gamma(i)=int_{t_i}^{t_{i+1}} (lambda_{i-1}B_{i-1}(t) +
 % lambda_iB_i(t))_+ B_i(t)dt,  i=2, ... N
 % Psi(i-1)=int_{t_i}^{t_{i+1}} (lambda_{i-1}B_{i-1}(t) +
 % lambda_iB_i(t))_+ B_{i-1}(t)dt,  i=2, ... N
 
 % for II(i)==-1 ([t_i, t_{i+1}] \in \Omgea_2
 % Gamma(i)= - int_{t_i}^{t_{i+1}} (lambda_{i-1}B_{i-1}(t) +
 % lambda_iB_i(t))_- B_i(t)dt,  i=2, ... N
 % Psi(i-1)= - int_{t_i}^{t_{i+1}} (lambda_{i-1}B_{i-1}(t) +
 % lambda_iB_i(t))_- B_{i-1}(t)dt,  i=2, ... N
 
 % for II(i)==0 ([t_i, t_{i+1}] \in \Omgea_3
 % Gamma(i)=int_{t_i}^{t_{i+1}} (lambda_{i-1}B_{i-1}(t) +
 % lambda_iB_i(t)) B_i(t)dt,  i=2, ... N
 % Psi(i-1)=int_{t_i}^{t_{i+1}} (lambda_{i-1}B_{i-1}(t) +
 % lambda_iB_i(t)) B_{i-1}(t)dt,  i=2, ... N
 
 Gamma=zeros(N,1);
 Psi=zeros(N,1);
 
 % Calculate rho(i):=B_i(t_{i+1}) and lrho(i):=lambda_i*rho(i)
 rho=zeros(N,1);
 lrho=zeros(N,1);
 for i=1:N
     rho(i)=2/(T(i+2)-T(i));
     lrho(i)=Lam(i)*rho(i);
 end
 
% Calculate the first piece of intergral
% [t_i, t_{i+1}]
%
  Phi1=0;
  if (II(1)==1 & Lam(1)>=0) | (II(1)==-1 & Lam(1)<=0) | (II(1)==0)
     h=(T(2)-T(1))/2;
     f=f+2*lrho(1)^2*h/3;
     Phi1=2*Lam(1)*rho(1)^2*h/3; % To calculate Phi_1
  end
  
% Calculate the pieces on [T(2), T(3)] ... [T(N), T(N+1)]
%
  for i=2:N
      if II(i)==1
         if (Lam(i-1)>=0) & (Lam(i) >=0)
            h=(T(i+1)-T(i))/2;
            f=f+( lrho(i-1)^2 + lrho(i)^2 + (lrho(i-1)+lrho(i))^2 )*h/3;
            Gamma(i)=(2*Lam(i)*rho(i)^2 + Lam(i-1)*rho(i-1)*rho(i))*h/3;
            Psi(i-1)=(2*Lam(i-1)*rho(i-1)^2+Lam(i)*rho(i-1)*rho(i))*h/3;
         else
            T0=(lrho(i)*T(i)-lrho(i-1)*T(i+1))/(lrho(i)-lrho(i-1));
            H=T(i+1)-T(i);
               if (Lam(i-1)<0) & (Lam(i)>0)
                  h=(T(i+1)-T0)/2;
                  f=f+2*lrho(i)^2*h/3;
                  b1=rho(i-1)*h/H; 
                  % b1=B_{i-1}((T(i+1)+T0)/2) middle poin value of B_{i-1} at
                  % (T(i+1)+T0)/2
                  b2=0.5*rho(i)*(1+(T0-T(i))/H);
                  % b2=B_i((T(i+1)+T0)/2) middle poin value of B_i at
                  % (T(i+1)+T0)/2
                  Gamma(i)=( Lam(i)*rho(i)^2 + 4*( Lam(i-1)*b1*b2 + Lam(i)*b2^2 ) )*h/3;
                  Psi(i-1)=4*( Lam(i-1)*b1^2 + Lam(i)*b1*b2 ) *h/3;
               end
               if (Lam(i-1)>0) & (Lam(i)<0)
                   h=(T0-T(i))/2;
                   f=f+2*lrho(i-1)^2*h/3;
                   b1=0.5*rho(i-1)*(1+(T(i+1)-T0)/H);
                    % b1=B_{i-1}((T(i)+T0)/2) middle poin value of B_{i-1} at
                    % (T(i)+T0)/2
                   b2=rho(i)*h/H;
                    % b2=B_i((T(i+1)+T0)/2) middle poin value of B_i at
                    % (T(i)+T0)/2
                   Gamma(i)=4*( Lam(i-1)*b1*b2 + Lam(i)*b2^2 )*h/3;
                   Psi(i-1)=( Lam(i-1)*rho(i-1)^2 + 4*(Lam(i-1)*b1^2+Lam(i)*b1*b2) )*h/3;
               end
         end 
      end % end of (II(i)==1)
         if II(i)==-1
            if (Lam(i-1)<=0) & (Lam(i) <=0)
            h=(T(i+1)-T(i))/2;
            f=f+( lrho(i-1)^2 + lrho(i)^2 + (lrho(i-1)+lrho(i))^2 )*h/3;
            Gamma(i)=(2*Lam(i)*rho(i)^2 + Lam(i-1)*rho(i-1)*rho(i))*h/3;
            Psi(i-1)=(2*Lam(i-1)*rho(i-1)^2+Lam(i)*rho(i-1)*rho(i))*h/3;
         else
            T0=(lrho(i)*T(i)-lrho(i-1)*T(i+1))/(lrho(i)-lrho(i-1));
            H=T(i+1)-T(i);
               if (Lam(i-1)>0) & (Lam(i)<0)
                  h=(T(i+1)-T0)/2;
                  f=f+2*lrho(i)^2*h/3;
                  b1=rho(i-1)*h/H; 
                  % b1=B_{i-1}((T(i+1)+T0)/2) middle poin value of B_{i-1} at
                  % (T(i+1)+T0)/2
                  b2=0.5*rho(i)*(1+(T0-T(i))/H);
                  % b2=B_i((T(i+1)+T0)/2) middle poin value of B_i at
                  % (T(i+1)+T0)/2
                  Gamma(i)=( Lam(i)*rho(i)^2 + 4*( Lam(i-1)*b1*b2 + Lam(i)*b2^2 ) )*h/3;
                  Psi(i-1)=4*( Lam(i-1)*b1^2 + Lam(i)*b1*b2 ) *h/3;
               end
               if (Lam(i-1)<0) & (Lam(i)>0)
                   h=(T0-T(i))/2;
                   f=f+2*lrho(i-1)^2*h/3;
                   b1=0.5*rho(i-1)*(1+(T(i+1)-T0)/H);
                    % b1=B_{i-1}((T(i)+T0)/2) middle poin value of B_{i-1} at
                    % (T(i)+T0)/2
                   b2=rho(i)*h/H;
                    % b2=B_i((T(i+1)+T0)/2) middle poin value of B_i at
                    % (T(i)+T0)/2
                   Gamma(i)=4*( Lam(i-1)*b1*b2 + Lam(i)*b2^2 )*h/3;
                   Psi(i-1)=( Lam(i-1)*rho(i-1)^2 + 4*(Lam(i-1)*b1^2+Lam(i)*b1*b2) )*h/3;
               end  
            end
         end % end of (II(i)==1)
         if II(i)==0
            h=(T(i+1)-T(i))/2;
            f=f+( lrho(i-1)^2 + lrho(i)^2 + (lrho(i-1)+lrho(i))^2 )*h/3;
            Gamma(i)=(2*Lam(i)*rho(i)^2 + Lam(i-1)*rho(i-1)*rho(i))*h/3;
            Psi(i-1)=(2*Lam(i-1)*rho(i-1)^2+Lam(i)*rho(i-1)*rho(i))*h/3;
         end % end of (II(i)==0)
  end % end for i=2:N

% Calculate the last piece of the integral on [t_{N+1}, t_{N+2}]
%
  Phi2=0;
  if (II(N)==1 & Lam(N)>=0) | (II(N)==0 & Lam(N)<=0) | (II(N)==0)
     h=(T(N+2)-T(N+1))/2;
     f=f+2*lrho(N)^2*h/3;
     Phi2=2*Lam(N)*rho(N)^2*h/3; % to calculate Phi_2
  end
% Calculate f
f=0.5*f - Lam'*D;

% Calculate F
 F=zeros(N,1);
 F(1)=Phi1+Psi(1);
 for i=2:(N-1)
     F(i)=Gamma(i)+Psi(i);
 end
 F(N)=Gamma(N)+Phi2;
 F=F-D;
 
 return;
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 %% Enf of Ff.m %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
  function [d, offd]=Jacob(T, Lam, II)

%The function fvalue uses one application of SIMPSON's rule to calculate
%numerically the integral which reaults in the Jacobian $M(Lam)$
%
% M(\lambda)(i,j) =\int_a^b (-\sum_{l=1}^N)_-^0 B_iB_j dt
% Refer to formulae (25) in my paper
% Output
% d=diagonal vector of the Jacobian matrix M
% offd=off diagonal vector (size (N-1)) of M

 [N,NC]=size(T);
 N=N-2;

% Calculate rho(i):=B_i(t_{i+1}) and lrho(i):=lambda_i*rho(i)
 rho=zeros(N,1);
 lrho=zeros(N,1);
 for i=1:N
     rho(i)=2/(T(i+2)-T(i));
     lrho(i)=Lam(i)*rho(i);
 end
 
 tau_1=0; tau_2=0;
 if (II(1)==1 & Lam(1)>=0) | (II(1)==-1 & Lam(1)<=0) | (II(1)==0)
     h=(T(2)-T(1))/2;
     tau_1=2*rho(1)^2*h/3;
 end
 
 if (II(N)==1 & Lam(N)>=0) | (II(N)==0 & Lam(N)<=0) | (II(N)==0)
     h=(T(N+2)-T(N+1))/2;
     tau_2=2*rho(N)^2*h/3;
 end
 % use a, b and c to denote the vectors contaning
 % a=(a_i) 
 % a_i=int_{t_i}^{t_{i+1}} P(lambda, t) B_i^2 dt
 % b=(b_i)
 % b_i=int_{t_i}^{t_{i+1}} P(lambda, t) B_{i-1}^2 dt
 % c=(c_i)
 % c_i=int_{t_i}^{t_{i+1}} P(lambda, t) B_iB_{i-1} dt
 % i=2, \ldots, N
 
 a=zeros(N-1,1);
 b=a;
 c=a;
 
 for i=2:N
     if II(i)==1
        if (Lam(i-1)>=0) & (Lam(i) >=0)
         h=(T(i+1)-T(i))/2;
         a(i-1)=2*rho(i)^2*h/3;
         b(i-1)=2*rho(i-1)^2*h/3;
         c(i-1)=rho(i-1)*rho(i)*h/3;
         else
         T0=(lrho(i)*T(i)-lrho(i-1)*T(i+1))/(lrho(i)-lrho(i-1));
         H=T(i+1)-T(i);
         biT0=rho(i)*(T0-T(i))/H;
         % biT0=B_i(T0)
         bi_1T0=rho(i-1)*(T(i+1)-T0)/H;
         % bi_1T0=B_{i-1}(T0)
         if (Lam(i-1)<0) & (Lam(i)>0)
             h=(T(i+1)-T0)/2;
             
             b1=rho(i-1)*h/H; 
             % b1=B_{i-1}((T(i+1)+T0)/2) middle poin value of B_{i-1} at
             % (T(i+1)+T0)/2
             b2=0.5*rho(i)*(1+(T0-T(i))/H);
             % b2=B_i((T(i+1)+T0)/2) middle poin value of B_i at
             % (T(i+1)+T0)/2
             a(i-1)=(rho(i)^2+biT0^2+4*b2^2)*h/3;
             b(i-1)=(bi_1T0^2+4*b1^2)*h/3;
             c(i-1)=(bi_1T0*biT0+4*b1*b2)*h/3;
         end
         if (Lam(i-1)>0) & (Lam(i)<0)
             h=(T0-T(i))/2;
             
             b1=0.5*rho(i-1)*(1+(T(i+1)-T0)/H);
             % b1=B_{i-1}((T(i)+T0)/2) middle poin value of B_{i-1} at
             % (T(i)+T0)/2
             b2=rho(i)*h/H;
             % b2=B_i((T(i)+T0)/2) middle poin value of B_i at
             % (T(i)+T0)/2
             a(i-1)=(biT0^2+4*b2^2)*h/3;
             b(i-1)=(rho(i-1)^2+4*b1^2+bi_1T0^2)*h/3;
             c(i-1)=(4*b1*b2+bi_1T0*biT0)*h/3;
         end
        end
     end % if II(i)==1
     if II(i)==-1
        if (Lam(i-1)<=0) & (Lam(i) <=0)
         h=(T(i+1)-T(i))/2;
         a(i-1)=2*rho(i)^2*h/3;
         b(i-1)=2*rho(i-1)^2*h/3;
         c(i-1)=rho(i-1)*rho(i)*h/3;
         else
         T0=(lrho(i)*T(i)-lrho(i-1)*T(i+1))/(lrho(i)-lrho(i-1));
         H=T(i+1)-T(i);
         biT0=rho(i)*(T0-T(i))/H;
         % biT0=B_i(T0)
         bi_1T0=rho(i-1)*(T(i+1)-T0)/H;
         % bi_1T0=B_{i-1}(T0)
         if (Lam(i-1)>0) & (Lam(i)<0)
             h=(T(i+1)-T0)/2;
             
             b1=rho(i-1)*h/H; 
             % b1=B_{i-1}((T(i+1)+T0)/2) middle poin value of B_{i-1} at
             % (T(i+1)+T0)/2
             b2=0.5*rho(i)*(1+(T0-T(i))/H);
             % b2=B_i((T(i+1)+T0)/2) middle poin value of B_i at
             % (T(i+1)+T0)/2
             a(i-1)=(rho(i)^2+biT0^2+4*b2^2)*h/3;
             b(i-1)=(bi_1T0^2+4*b1^2)*h/3;
             c(i-1)=(bi_1T0*biT0+4*b1*b2)*h/3;
         end
         if (Lam(i-1)<0) & (Lam(i)>0)
             h=(T0-T(i))/2;
             
             b1=0.5*rho(i-1)*(1+(T(i+1)-T0)/H);
             % b1=B_{i-1}((T(i)+T0)/2) middle poin value of B_{i-1} at
             % (T(i)+T0)/2
             b2=rho(i)*h/H;
             % b2=B_i((T(i)+T0)/2) middle poin value of B_i at
             % (T(i)+T0)/2
             a(i-1)=(biT0^2+4*b2^2)*h/3;
             b(i-1)=(rho(i-1)^2+4*b1^2+bi_1T0^2)*h/3;
             c(i-1)=(4*b1*b2+bi_1T0*biT0)*h/3;
         end
        end 
     end % if II(i)==-1
     if II(i)==0
        h=(T(i+1)-T(i))/2;
        a(i-1)=2*rho(i)^2*h/3;
        b(i-1)=2*rho(i-1)^2*h/3;
        c(i-1)=rho(i-1)*rho(i)*h/3; 
     end % if II(i)==0     
 end % for i=2:N

% Assigin values to d and offd
  d=zeros(N,1);
  offd=zeros(N-1,1);
  
  d(1)=tau_1+b(1);
  offd(1)=c(1);
  d(N)=tau_2+a(N-1);
  for i=2:(N-1)
      d(i)=a(i-1)+b(i);
      offd(i)=c(i);
  end
 
  return;
  
%%%%%%%%%%%%%% Enf of Jacob.m
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function b = sptsl ( n, d, e, b )

%% SPTSL solves a positive definite tridiagonal linear system.
%
%  Modified:
%
%    09 November 2006
%
%  Author:
%
%    MATLAB translation by John Burkardt.
%
%  Reference:
%
%    Dongarra, Moler, Bunch and Stewart,
%    LINPACK User's Guide,
%    SIAM, (Society for Industrial and Applied Mathematics),
%    3600 University City Science Center,
%    Philadelphia, PA, 19104-2688.
%    ISBN 0-89871-172-X
%
%  Parameters:
%
%    Input, integer N, the order of the matrix.
%
%    Input, real D(N), the diagonal of the tridiagonal matrix.
%
%    Input, E(N), the offdiagonal of the tridiagonal matrix in
%    entries E(1:N-1).
%
%    Input, real B(N), the right hand side.
%
%    Output, real B(N), the solution.
%

%
%  Check for 1 x 1 case.
%
  if ( n == 1 )
    b(1) = b(1) / d(1);
    return
  end

  nm1d2 = floor ( ( n - 1 ) / 2 );

  if ( 2 < n )

    kbm1 = n - 1;
%
%  Zero top half of subdiagonal and bottom half of superdiagonal.
%
    for k = 1 : nm1d2
      t1 = e(k) / d(k);
      d(k+1) = d(k+1) - t1 * e(k);
      b(k+1) = b(k+1) - t1 * b(k);
      t2 = e(kbm1) / d(kbm1+1);
      d(kbm1) = d(kbm1) - t2 * e(kbm1);
      b(kbm1) = b(kbm1) - t2 * b(kbm1+1);
      kbm1 = kbm1 - 1;
    end

  end

  kp1 = nm1d2 + 1;
%
%  Clean up for possible 2 x 2 block at center.
%
  if ( mod ( n, 2 ) == 0 )
    t1 = e(kp1) / d(kp1);
    d(kp1+1) = d(kp1+1) - t1 * e(kp1);
    b(kp1+1) = b(kp1+1) - t1 * b(kp1);
    kp1 = kp1 + 1;
  end
%
%  Back solve starting at the center, going towards the top and bottom.
%
  b(kp1) = b(kp1) / d(kp1);

  if ( 2 < n )

    k = kp1 - 1;
    ke = kp1 + nm1d2 - 1;

    for kf = kp1 : ke
      b(k) = ( b(k) - e(k) * b(k+1) ) / d(k);
      b(kf+1) = ( b(kf+1) - e(kf) * b(kf) ) / d(kf+1);
      k = k - 1;
    end

  end

  if ( mod ( n, 2 ) == 0 )
    b(1) = ( b(1) - e(1) * b(2) ) / d(1);
  end

return;

%%%%% End of sptsl.m %%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function cubic=Cubic(T,Y,II,alpha)
% II is the vector indicating where the intervals begong
% (T, Y) are original data
% alpha is the solution vector of the Dnewton
% alpha=Lam

  [N, Nc]=size(T);
  N=N-2;


%  for i=1:N
%    alpha(i)=alpha(i)/(T(i+2)-T(i));
%  end

% Between two endpoints, there exists at most one new Knot
% So for every interval [t_i, t_{i+1}], we assign six indexes, indicating
% where the knots are located
% K(i,j), j=1,2 3. and Z(i,j), j=1, 2,3.
% K(i,1)=t_i, K(i,3)=t_{i+1}, K(i,2) be the possible new Knot
% Z(i,1), Z(i,2) and Z(i,3) be the second derivatives at the Knots
% K(i,j), j=1,2,3.

% Sampling on each interval [T_i, T_{i+1}] and calculating the functional
% values on these sampling points. In every interval, we take 10 sample
% points, i.e., the jump step is h(i,1)/10 or h(i,2)/10 below.
%
% xs is x-sampling vector while ys is y-sampling vector

  xs=[]; ys=[];
  
% Working Interval [t_1, t_2]

  if II(1)==1
     K(1,1)=T(1); K(1,2)=0; K(1,3)=T(2);
     Z(1,1)=0;    Z(1,2)=0; Z(1,3)=max(0,alpha(1));
  elseif II(1)==-1
     K(1,1)=T(1); K(1,2)=0; K(1,3)=T(2);
     Z(1,1)=0;    Z(1,2)=0; Z(1,3)=min(0,alpha(1));
  end

% Cubic Spline in the interval [t_1, t_2]
% S(1,2)(t)=Z(1,3)/(6*h(1,2))*(t-t(1))^3+C(1,2)(t-K(1,1))+D(1,2)(K(1,3)-t)
% Conditions: S(1,2)(T(1))=Y(1);  S(1,2)(T_2)=Y(2)

  C(1,1)=0; D(1,1)=0;   % No use
  h(1,1)=0; %no use
  h(1,2)=K(1,3)-K(1,1);  %=T(2)-T(1);
  C(1,2)=Y(2)/h(1,2)-Z(1,3)*h(1,2)/6;
  D(1,2)=Y(1)/h(1,2);

% Save the first derivative of S(1,2) at t_2 for later use, denoted by ds

  ds=Z(1,3)/2*h(1,2)+C(1,2)-D(1,2);
 
% Sampling

  xsample=K(1,1):h(1,2)/10:K(1,3);
  ysample=Z(1,3)/(6*h(1,2))*(xsample-K(1,1)).^3+C(1,2)*(xsample-K(1,1))+...
          D(1,2)*(K(1,3)-xsample);
  xs=[xs xsample]; ys=[ys ysample];
  
% Now working on intervals [t_2, t_3] .... [t_N, t_{N+1}]

  for i=2:N

    if II(i)==1

% if min(alpha(i-1), alpha(i))<0 & max(alpha(i-1), alpha(i))>0    
% There is a new knot inserted
% Otherwise, no new knot incresed

  if alpha(i-1)>0 & alpha(i)<0
     K(i,1)=T(i); 
     K(i,2)=(alpha(i)*T(i)-alpha(i-1)*T(i+1))/(alpha(i)-alpha(i-1));
     K(i,3)=T(i+1);

     Z(i,1)=alpha(i-1);
     Z(i,2)=0;
     Z(i,3)=0;
 
% Cubic spline: There are two cubic splines on the interval 
% [T_i, K(i,2)] and [K(i,2), T_{i+1}], respectively
% Working interval [T_i, K(i,2)]
% S(i,1)(t)= Z(i,1)/(6*h(i,1))*(K(i,2)-t)^3+C(i,1)(t-K(i,1))+D(i,1)(K(i,2)-t)
% h(i,1)=K(i,2)-K(i,1)
% Conditions S(i,1)(T_i)=Y(i), S(i,1)'(T_i)=S(i-1,2)'(T_i)

  h(i,1)=K(i,2)-K(i,1);
  D(i,1)=Y(i)/h(i,1)-Z(i,1)*h(i,1)/6;
  C(i,1)=ds+Z(i,1)*h(i,1)/2+D(i,1);

% Save the function value at K(i,2)=yy;

  yy=C(i,1)*h(i,1);

% Sampling

  xsample=K(i,1):h(i,1)/10:K(i,2);
  ysample=Z(i,1)/(6*h(i,1))*(K(i,2)-xsample).^3+C(i,1)*(xsample-K(i,1))+...
  D(i,1)*(K(i,2)-xsample);
  xs=[xs xsample]; ys=[ys ysample];

% Working Interval [K(i,2), K(i,3)]
% S(i,2)(t)=C(i,2)(t-K(i,2))+D(i,2)(K(i,3)-t)
% Conditions: S(i,2)(T_{i+1})=Y(i+1); S(i,2)(K(i,2)=y

  h(i,2)=T(i+1)-K(i,2);
  C(i,2)=Y(i+1)/h(i,2);
  D(i,2)=yy/h(i,2);

% % Save the first derivative of S(i,2) at t_{i+1} for later use, denoted by ds

  ds=C(i,2)-D(i,2);

% Sampling

  xsample=K(i,2):h(i,2)/10:K(i,3);
  ysample=C(i,2)*(xsample-K(i,2))+D(i,2)*(K(i,3)-xsample);
  xs=[xs xsample]; ys=[ys ysample];

  elseif alpha(i-1)<0 & alpha(i)>0
  
    K(i,1)=T(i);    Z(i,1)=0;
    K(i,2)=(alpha(i)*T(i)-alpha(i-1)*T(i+1))/(alpha(i)-alpha(i-1));
    Z(i,2)=0;
    K(i,3)=T(i+1);  Z(i,3)=alpha(i);
    
% Cubic Splines: There are two cubic splines on the inetrvals:
% [K(i,1) K(i,2)] and [K(i,2), K(i,3)], respectively
% Working interval [K(i,1), K(i,2)]
% S(i,1)(t)=C(i,1)(t-K(i,1))+D(i,1)(K(i,2)-t)
% h(i,1)=K(i,2)-K(i,1);
% Conditions: S(i,1)(T(i))=Y(i), S(i,1)'(T(i))=S(i-1,2)'(T(i))

  h(i,1)=K(i,2)-K(i,1);
  D(i,1)=Y(i)/h(i,1);
  C(i,1)=ds+D(i,1);
  
% Save the function value at K(i,2), denoted by yy

  yy=C(i,1)*h(i,1);
 
% Sampling

  xsample=K(i,1):h(i,1)/10:K(i,2);
  ysample=C(i,1)*(xsample-K(i,1))+D(i,1)*(K(i,2)-xsample);
  xs=[xs xsample]; ys=[ys ysample];
  
% Working Interval [K(i,2) K(i,3)]
% S(i,2)(t)=Z(i,3)/(6*h(i,2))*(t-K(i,2))^3+C(i,2)*(t-K(i,2))+D(i,2)*(K(i,3)-t)
% h(i,2)=K(i,3)-K(i,2);
% Conditions: S(i,2)(K(i,2))=yy, S(i,2)(T(i+1))=Y(i+1)

  h(i,2)=K(i,3)-K(i,2);
  D(i,2)=yy/h(i,2);
  C(i,2)=Y(i+1)/h(i,2)-Z(i,3)*h(i,2)/6;
  
% Save the first derivative of S(i,2) at t_{i+1} for later use, denoted by ds 
 
  ds=Z(i,3)*h(i,2)/2+C(i,2)-D(i,2);
  
% Sampling
  
  xsample=K(i,2):h(i,2)/10:K(i,3);
  ysample=Z(i,3)/(6*h(i,2))*(xsample-K(i,2)).^3+...
          C(i,2)*(xsample-K(i,2))+D(i,2)*(K(i,3)-xsample);
  xs=[xs xsample]; ys=[ys ysample];
   
  else   % No new knots added 

    K(i,1)=T(i);   Z(i,1)=max(0, alpha(i-1));
    K(i,2)=0;      Z(i,2)=0;  %No use
    K(i,3)=T(i+1); Z(i,3)=max(0,alpha(i));

% Cubic Spline
% S(i,2)(t)=Z(i,3)/(6h(i,2))(t-K(i,1))^3+Z(i,1)/(6h(i,2))(K(i,3)-t)^3+...
%           C(i,2)(t-K(i,1))+D(i,2)(K(i,3)-t);
% Only one cubic spline, denoted S(i,2)(t)
% Conditions: S(i,2)(T_i)=Y_i; S(i,2)(T_{i+1})=Y_{i+1}

   h(i,2)=T(i+1)-T(i);
   D(i,2)=Y(i)/h(i,2)-Z(i,1)/6*h(i,2);
   C(i,2)=Y(i+1)/h(i,2)-Z(i,3)/6*h(i,2);
% Save the first derivative of S(i,2) at t_{i+1} for later use, denoted by ds

  ds=Z(i,3)/2*h(i,2)+C(i,2)-D(i,2);

% Sampling

  xsample=K(i,1):h(i,2)/10:K(i,3);
  ysample=Z(i,3)/(6*h(i,2))*(xsample-K(i,1)).^3+...
          Z(i,1)/(6*h(i,2))*(K(i,3)-xsample).^3+...
          C(i,2)*(xsample-K(i,1))+D(i,2)*(K(i,3)-xsample);
  xs=[xs xsample]; ys=[ys ysample];

  end
 end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%

  if II(i)==-1

% if min(alpha(i-1), alpha(i))<0 & max(alpha(i-1), alpha(i))>0
% There is a new knot inserted
% Otherwise, no new knot incresed
 
     if alpha(i-1)>0 & alpha(i)<0

        K(i,1)=T(i);   Z(i,1)=0;
        K(i,2)=(alpha(i)*T(i)-alpha(i-1)*T(i+1))/(alpha(i)-alpha(i-1));
        Z(i,2)=0;
        K(i,3)=T(i+1); Z(i,3)=alpha(i);

% Cubic spline: There are two cubic splines on [T(i), K(i,2)] and 
% [K(i,2), T(i+1)], respectively
% Working Interval [T(i), K(i,2)]
% S(i,1)(t)=C(i,1)(t-K(i,1))+D(i,1)(K(i,2)-t);
% Conditions: S(i,1)(T(i))=Y(i), S(i,1)'(T(i))=S(i-1,2)'(T(i));

        h(i,1)=K(i,2)-T(i);
        D(i,1)=Y(i)/h(i,1);
        C(i,1)=ds+D(i,1);

% Save the function value at K(i,2)=yy;
        
        yy=C(i,1)*h(i,1);

% Sampling

  xsample=K(i,1):h(i,1)/10:K(i,2);
  ysample=C(i,1)*(xsample-K(i,1))+D(i,1)*(K(i,2)-xsample);
  xs=[xs xsample]; ys=[ys ysample];

% Working Interval [K(i,2), T(i+1)]=[K(i,2), K(i,3)] 
% S(i,2)(t)=Z(i,3)/(6h(i,2))(t-K(i,2))^3+C(i,2)(t-K(i,2))+D(i,2)(K(i,3)-t);
%Conditions: S(i,2)(K(i,2))=y; S(i,2)(K(i,3))=Y(i+1);

        h(i,2)=K(i,3)-K(i,2);
        D(i,2)=yy/h(i,2);
        C(i,2)=Y(i+1)/h(i,2)-Z(i,3)/6*h(i,2);
% Save the first derivative of S(i,2) at t_{i+1} for later use, denoted by ds

        ds=Z(i,3)/2*h(i,2)+C(i,2)-D(i,2);
% Sampling

  xsample=K(i,2):h(i,2)/10:K(i,3);
  ysample=Z(i,3)/(6*h(i,2))*(xsample-K(i,2)).^3+...
          C(i,2)*(xsample-K(i,2))+D(i,2)*(K(i,3)-xsample);
  xs=[xs xsample]; ys=[ys ysample];

      elseif alpha(i-1)<0 & alpha(i)>0
      
      K(i,1)=T(i);     Z(i,1)=alpha(i-1);
      K(i,2)=(alpha(i)*T(i)-alpha(i-1)*T(i+1))/(alpha(i)-alpha(i-1));
      Z(i,2)=0;
      K(i,3)=T(i+1);   Z(i,3)=0;
      
% Cubic spline: There are two cubic splines on 
% [K(i,1), K(i,2)] and [K(i,2), K(i,3)], respectively
% Working Interval [K(i,1), K(i,2)]
% S(i,1)=Z(i,1)/(6h(i,1))*(K(i,2)-t)^3+C(i,1)(t-K(i,1))+D(i,1)(K(i,2)-t)
% h(i,1)=K(i,2)-K(i,1)
% Conditions: S(i,1)(T(i))=Y(i), S(i,1)'(T(i))=S(i-1,2)'(T(i))=ds

      h(i,1)=K(i,2)-K(i,1);
      D(i,1)=Y(i)/h(i,1)-Z(i,1)*h(i,1)/6;
      C(i,1)=ds+D(i,1)+Z(i,1)*h(i,1)/2;
      
% Save the function value at $K(i,2), denoted by yy

      yy=C(i,1)*h(i,1);
      
% Sampling 
    
      xsample=K(i,1):h(i,1)/10:K(i,2);
      ysample=Z(i,1)/(6*h(i,1))*(K(i,2)-xsample).^3+...
              C(i,1)*(xsample-K(i,1))+D(i,1)*(K(i,2)-xsample);
      xs=[xs xsample]; ys=[ys ysample];
      
% Working interval" [K(i,2), K(i,3)]
% S(i,2)(t)=C(i,2)(t-K(i,2))+D(i,2)(K(i,3)-t);
% Conditiond: S(i,2)(K(i,2))=yy; S(i,2)(T(i+1))=Y(i+1);

      h(i,2)=K(i,3)-K(i,2);
      D(i,2)=yy/h(i,2);
      C(i,2)=Y(i+1)/h(i,2);
      
% Save the first derivative of S(i,2) at t_{i+1} for later use, denoted by ds
  
      ds=C(i,2)-D(i,2);
      
% Sampling

     xsample=K(i,2):h(i,2)/10:K(i,3);    
     ysample=C(i,2)*(xsample-K(i,2))+D(i,2)*(K(i,3)-xsample);
     xs=[xs xsample]; ys=[ys ysample];    
            
      else  % No New knots added

        K(i,1)=T(i);   Z(i,1)=min(0, alpha(i-1));
        K(i,2)=0;      Z(i,2)=0; %No use
        K(i,3)=T(i+1); Z(i,3)=min(0,alpha(i));

% There is only one piece of cubic spline, namely
% S(i,2)(t)=Z(i,3)/(6h(i,2))(t-K(i,1))^3+Z(i,1)/(6h(i,2))(K(i,3)-t)^3+...
%           C(i,2)(t-K(i,1))+D(i,2)(K(i,3)-t);
% Conditions: S(i,2)(T(i))=Y(i);  S(i,2)(T(i+1))=Y(i+1);

        h(i,2)=K(i,3)-K(i,1);
        D(i,2)=Y(i)/h(i,2)-Z(i,1)/6*h(i,2);
        C(i,2)=Y(i+1)/h(i,2)-Z(i,3)/6*h(i,2);
% Save the first derivative of S(i,2) at t_{i+1} for later use, denoted by ds
     
        ds=Z(i,3)*h(i,2)/2+C(i,2)-D(i,2); 
% Sampling

  xsample=K(i,1):h(i,2)/10:K(i,3);
  ysample=Z(i,3)/(6*h(i,2))*(xsample-K(i,1)).^3+...
          Z(i,1)/(6*h(i,2))*(K(i,3)-xsample).^3+...
          C(i,2)*(xsample-K(i,1))+D(i,2)*(K(i,3)-xsample);
  xs=[xs xsample]; ys=[ys ysample];
     end
  end
    
    if II(i)==0
%This is natural spline on [T_i, T_{i+1}]

        K(i,1)=T(i);   Z(i,1)=alpha(i-1);
        K(i,2)=0;      Z(i,2)=0;   %No use
        K(i,3)=T(i+1); Z(i,3)=alpha(i);       
% Cubic spline: Only one, denoted by S(i,2)(t)
% S(i,2)(t)=Z(i,3)/(6h(i,2))(t-K(i,1))^3+Z(i,1)/(6h(i,2))(K(i,3)-t)^3+...
%           C(i,2)(t-K(i,1))+D(i,2)(K(i,3)-t);
%Conditions: S(i,2)(T(i))=Y(i);  S(i,2)(T(i+1))=Y(i+1)

        h(i,2)=T(i+1)-T(i);
        D(i,2)=Y(i)/h(i,2)-Z(i,1)/6*h(i,2);
        C(i,2)=Y(i+1)/h(i,2)-Z(i,3)/6*h(i,2);
% Save the first derivative of S(i,2) at t_{i+1} for later use, denoted by ds
     
        ds=Z(i,3)*h(i,2)/2+C(i,2)-D(i,2);
% Sampling

  xsample=K(i,1):h(i,2)/10:K(i,3);
  ysample=Z(i,3)/(6*h(i,2))*(xsample-K(i,1)).^3+...
          Z(i,1)/(6*h(i,2))*(K(i,3)-xsample).^3+...
          C(i,2)*(xsample-K(i,1))+D(i,2)*(K(i,3)-xsample);
  xs=[xs xsample]; ys=[ys ysample];
  end
end
%Working on the last interval [T(N+1), T(N+2)]
% This interval is determined by II(N), i.e., II(N)==1, it is convex ...

    if II(N)==1
       K(N+1,1)=T(N+1);   Z(N+1,1)=max(0,alpha(N));
       K(N+1,2)=0;        Z(N+1,2)=0;  %No use
       K(N+1,3)=T(N+2);   Z(N+1,3)=0;
     elseif II(N)==-1
       K(N+1,1)=T(N+1);   Z(N+1,1)=min(0,alpha(N));
       K(N+1,2)=0;        Z(N+1,2)=0;  %No use
       K(N+1,3)=T(N+2);   Z(N+1,3)=0;
     else
       K(N+1,1)=T(N+1);   Z(N+1,1)=alpha(N);
       K(N+1,2)=0;        Z(N+1,2)=0;  %No use
       K(N+1,3)=T(N+2);   Z(N+1,3)=0; 
    end
%Only one piece of cubic spline
%Denoted S(N+1,2)(t)=Z(N+1,1)/(6h(N+1,2))(T(N+2)-t)^3+C(N+1,2)(t-T(N+1))+...
%                    D(N+1,2)(T(N+2)-t);
%Conditions: S(N+1,2)(T(N+1))=Y(N+1);  S(N+1,2)(T(N+2))=Y(N+2)

    h(N+1,2)=T(N+2)-T(N+1);
    D(N+1,2)=Y(N+1)/h(N+1,2)-Z(N+1,1)*h(N+1,2)/6;
    C(N+1,2)=Y(N+2)/h(N+1,2);

% Sampling

  xsample=K(N+1,1):h(N+1,2)/10:K(N+1,3);
  ysample=Z(N+1,1)/(6*h(N+1,2))*(T(N+2)-xsample).^3+...
          C(N+1,2)*(xsample-T(N+1))+D(N+1,2)*(T(N+2)-xsample);
  xs=[xs xsample]; ys=[ys ysample]; 

% Generate the natural cubic spline
%
  xxx=T(1):0.02:T(N+2);
  yyy=spline(T,Y,xxx);
  
% Plot the natural cubic spline and the shape-preserving spline in 
% one figure

% plot(T,Y,'o',xxx,yyy,'-',xs,ys,'--')

% Plot the shape-preserving spline only 
plot(T,Y,'o',xs,ys,'--')
  
  cubic=1;

return;

%%%%%%%end of Cubic.m %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%