function [dbetaL] = Newton_Increment ( k_w1, k_w2, k_theta1, k_theta2, beta);

%  characteristic matrix for beam vibration problem;
%  the two end supports are constrained by two transverse spring
%  as well as two rotational springs;


%input beta = beta*L;

%%%%%%%%%%%%%%%%%%%
%  data preparation
%%%%%%%%%%%%%%%%%%%

cc = cos (beta);
ss = sin (beta);

ch = cosh(beta);
sh = sinh(beta);

Cbeam = [     
 beta^3                         -k_w1              -beta^3                                 -k_w1; 
 -k_theta1                     -beta               -k_theta1                                beta;			 
-(beta^3*cc+k_w2*ss)      (beta^3*ss-k_w2*cc)     (beta^3*ch-k_w2*sh)       (beta^3*sh-k_w2*ch);
 (beta*ss - k_theta2*cc)  (beta*cc + k_theta2*ss)  -(beta*sh +k_theta2*ch)   -(beta*ch +k_theta2*sh)];

% residual
residual = det(Cbeam); 

 C11 = [     
 3*beta^2                         -k_w1              -beta^3                                 -k_w1; 
 -k_theta1                     -beta               -k_theta1                                beta;			 
-(beta^3*cc+k_w2*ss)      (beta^3*ss-k_w2*cc)     (beta^3*ch-k_w2*sh)       (beta^3*sh-k_w2*ch);
 (beta*ss - k_theta2*cc)  (beta*cc + k_theta2*ss)  -(beta*sh +k_theta2*ch)   -(beta*ch +k_theta2*sh)];

 
 C12 = [     
 beta^3                         0               -beta^3                                 -k_w1; 
 -k_theta1                     -beta               -k_theta1                                beta;			 
-(beta^3*cc+k_w2*ss)      (beta^3*ss-k_w2*cc)     (beta^3*ch-k_w2*sh)       (beta^3*sh-k_w2*ch);
 (beta*ss - k_theta2*cc)  (beta*cc + k_theta2*ss)  -(beta*sh +k_theta2*ch)   -(beta*ch +k_theta2*sh)];


 C13 = [     
 beta^3                         -k_w1              -3*beta^2                                 -k_w1; 
 -k_theta1                     -beta               -k_theta1                                beta;			 
-(beta^3*cc+k_w2*ss)      (beta^3*ss-k_w2*cc)     (beta^3*ch-k_w2*sh)       (beta^3*sh-k_w2*ch);
 (beta*ss - k_theta2*cc)  (beta*cc + k_theta2*ss)  -(beta*sh +k_theta2*ch)   -(beta*ch +k_theta2*sh)];


 C14 = [     
 beta^3                         -k_w1              -beta^3                                    0; 
 -k_theta1                     -beta               -k_theta1                                beta;			 
-(beta^3*cc+k_w2*ss)      (beta^3*ss-k_w2*cc)     (beta^3*ch-k_w2*sh)       (beta^3*sh-k_w2*ch);
 (beta*ss - k_theta2*cc)  (beta*cc + k_theta2*ss)  -(beta*sh +k_theta2*ch)   -(beta*ch +k_theta2*sh)];

% newton slop
slope1 = det(C11) + det(C12) + det(C13) + det(C14);

 C21 = [     
 beta^3                         -k_w1              -beta^3                                 -k_w1; 
  0                            -beta               -k_theta1                                beta;			 
-(beta^3*cc+k_w2*ss)      (beta^3*ss-k_w2*cc)     (beta^3*ch-k_w2*sh)       (beta^3*sh-k_w2*ch);
 (beta*ss - k_theta2*cc)  (beta*cc + k_theta2*ss)  -(beta*sh +k_theta2*ch)   -(beta*ch +k_theta2*sh)];


 C22 = [     
 beta^3                         -k_w1              -beta^3                                 -k_w1; 
 -k_theta1                     -1.0               -k_theta1                                beta;			 
-(beta^3*cc+k_w2*ss)      (beta^3*ss-k_w2*cc)     (beta^3*ch-k_w2*sh)       (beta^3*sh-k_w2*ch);
 (beta*ss - k_theta2*cc)  (beta*cc + k_theta2*ss)  -(beta*sh +k_theta2*ch)   -(beta*ch +k_theta2*sh)];


 C23= [     
 beta^3                         -k_w1              -beta^3                                 -k_w1; 
 -k_theta1                     -beta                0                                      beta;			 
-(beta^3*cc+k_w2*ss)      (beta^3*ss-k_w2*cc)     (beta^3*ch-k_w2*sh)       (beta^3*sh-k_w2*ch);
 (beta*ss - k_theta2*cc)  (beta*cc + k_theta2*ss)  -(beta*sh +k_theta2*ch)   -(beta*ch +k_theta2*sh)];


 C24 = [     
 beta^3                         -k_w1              -beta^3                                 -k_w1; 
 -k_theta1                     -beta               -k_theta1                                1.0;			 
-(beta^3*cc+k_w2*ss)      (beta^3*ss-k_w2*cc)     (beta^3*ch-k_w2*sh)       (beta^3*sh-k_w2*ch);
 (beta*ss - k_theta2*cc)  (beta*cc + k_theta2*ss)  -(beta*sh +k_theta2*ch)   -(beta*ch +k_theta2*sh)];
% newton slope
slope2 = det(C21) + det(C22) + det(C23) + det(C24);


 C31 = [     
 beta^3                         -k_w1              -beta^3                                 -k_w1; 
 -k_theta1                     -beta               -k_theta1                                beta;			 
-(3*beta^2*cc-beta^3*ss +k_w2*cc)  (beta^3*ss-k_w2*cc)     (beta^3*ch-k_w2*sh)       (beta^3*sh-k_w2*ch);
 (beta*ss - k_theta2*cc)  (beta*cc + k_theta2*ss)  -(beta*sh +k_theta2*ch)   -(beta*ch +k_theta2*sh)];


 C32 = [     
 beta^3                         -k_w1              -beta^3                                 -k_w1; 
 -k_theta1                     -beta               -k_theta1                                beta;			 
-(beta^3*cc+k_w2*ss)      (3*beta^2*ss+ beta^3*cc +k_w2*ss)     (beta^3*ch-k_w2*sh)       (beta^3*sh-k_w2*ch);
 (beta*ss - k_theta2*cc)  (beta*cc + k_theta2*ss)  -(beta*sh +k_theta2*ch)   -(beta*ch +k_theta2*sh)];


 C33 = [     
 beta^3                         -k_w1              -beta^3                                 -k_w1; 
 -k_theta1                     -beta               -k_theta1                                beta;			 
-(beta^3*cc+k_w2*ss)      (beta^3*ss-k_w2*cc)     (3*beta^2*ch+beta^3*sh-k_w2*ch)       (beta^3*sh-k_w2*ch);
 (beta*ss - k_theta2*cc)  (beta*cc + k_theta2*ss)  -(beta*sh +k_theta2*ch)   -(beta*ch +k_theta2*sh)];

 
 C34 = [     
 beta^3                         -k_w1              -beta^3                                 -k_w1; 
 -k_theta1                     -beta               -k_theta1                                beta;			 
-(beta^3*cc+k_w2*ss)      (beta^3*ss-k_w2*cc)     (beta^3*ch-k_w2*sh)       (3*beta^2*sh+beta^3*ch-k_w2*sh);
 (beta*ss - k_theta2*cc)  (beta*cc + k_theta2*ss)  -(beta*sh +k_theta2*ch)   -(beta*ch +k_theta2*sh)];
% newton slope
slope3 = det(C31) + det(C32) + det(C33) + det(C34);

  
 C41 = [     
 beta^3                         -k_w1              -beta^3                                 -k_w1; 
 -k_theta1                     -beta               -k_theta1                                beta;			 
-(beta^3*cc+k_w2*ss)      (beta^3*ss-k_w2*cc)     (beta^3*ch-k_w2*sh)       (beta^3*sh-k_w2*ch);
 (ss + beta*cc + k_theta2*ss)  (beta*cc + k_theta2*ss)  -(beta*sh +k_theta2*ch)   -(beta*ch +k_theta2*sh)];

 
 C42 = [     
 beta^3                         -k_w1              -beta^3                                 -k_w1; 
 -k_theta1                     -beta               -k_theta1                                beta;			 
-(beta^3*cc+k_w2*ss)      (beta^3*ss-k_w2*cc)     (beta^3*ch-k_w2*sh)       (beta^3*sh-k_w2*ch);
 (beta*ss - k_theta2*cc)  (cc-beta*ss + k_theta2*cc)  -(beta*sh +k_theta2*ch)   -(beta*ch +k_theta2*sh)];

 
 C43 = [     
 beta^3                         -k_w1              -beta^3                                 -k_w1; 
 -k_theta1                     -beta               -k_theta1                                beta;			 
-(beta^3*cc+k_w2*ss)      (beta^3*ss-k_w2*cc)     (beta^3*ch-k_w2*sh)       (beta^3*sh-k_w2*ch);
 (beta*ss - k_theta2*cc)  (beta*cc + k_theta2*ss)  -(sh + beta*ch +k_theta2*sh)   -(beta*ch +k_theta2*sh)];

 
 C44 = [     
 beta^3                         -k_w1              -beta^3                                 -k_w1; 
 -k_theta1                     -beta               -k_theta1                                beta;			 
-(beta^3*cc+k_w2*ss)      (beta^3*ss-k_w2*cc)     (beta^3*ch-k_w2*sh)       (beta^3*sh-k_w2*ch);
 (beta*ss - k_theta2*cc)  (beta*cc + k_theta2*ss)  -(beta*sh +k_theta2*ch)   -(ch  + beta*sh +k_theta2*ch)];
% newton slope
slope4 = det(C41) + det(C42) + det(C43) + det(C44);


% compute Newton increment
slope = slope1 + slope2 + slope3  + slope4;
if abs(slope)<0.01
	dbetaL = 0;
else
    dbetaL = -	residual/slope;
end;	
return;