function Tss = ssnr(Tinit,fcorr,epsnr)
% Newton-Raphson iteration to find steady-state temperature profile 
% of reactor with preheater. (It may be unstable if we do not apply control)
%   Tinit - initial temperaure profile
%   fcorr - correction for NR-iteration; reduce if convergense problems
%   Tss = steady-state temperature profile
global mc mh Cp Cpc Q1 Q2 Q3 Npoint Tf cf dm1 dm2 dm3 dHrx u conc

% Linearize to find Jacobi matrix
delta=0.001;
for i=1:Npoint
 perturb=zeros(Npoint,1);
 perturb(i)=delta;
 pert1 = righthand(0,Tinit+perturb')';
 pert2 = righthand(0,Tinit-perturb')';
 A0(:,i)=( pert1 - pert2 )/(2*delta);
end
 J=inv(A0);

% Newton-Raphson iteration
 error=1;
 Tss=Tinit;
while error>epsnr
    func=righthand(0,Tss); 
    Tss=(Tss'-fcorr*J*func')';
    error=norm(func');
    error_steady_state_temperature_profile=error
  end

