N=100;                            %節點數
P0=linspace(0,0,N+1);            %初始壓力P0向量
P=linspace(0,0,N+1);             %壓力P向量
DP=linspace(0,0,N+1);            %ΔP
XP=linspace(0,0,N+1);            %各個節點值X
H=linspace(0,0,N+1);             %膜厚H
F=linspace(0,0,N+1);             %F向量
B=linspace(0,0,N+1);             %B向量
DF=linspace(0,0,N+1);            %ΔF
DL=linspace(0,0,N+1);            %求CIJ時的一個中間變量
MIDU=linspace(0,0,N+1);          %密度壓力關系
NIANDU=linspace(0,0,N+1);        %粘度壓力關系
CIJ=zeros(N+1,N+1);              %C矩陣
KK=zeros(N,N);               %最終的K矩陣N*N
KK1=zeros(N+1,N+1);              %最終的K1矩陣
KK2=zeros(N+1,N+1);              %最終的K2矩陣
LNODS=zeros(N+1,2);              %中間矩陣

initialization;         %初始化
for i=1:N+1
    if(XP(i)^2<=1)
        P0(i)=sqrt(1-XP(i)^2);
    end
end
P=P0;
DH0=1;
MAXIT=100;
ITER=0;
while (abs(DH0)>0.01&&ITER<=MAXIT)
    ITER=ITER+1;
    IT=0;
    EP=1;
    while (EP>0.01&&IT<=MAXIT)
        IT=IT+1;
        film_thick;          %求膜厚方程
        matrix_final;        %求解最終矩陣
       DP(2:N+1)=KK\DF(2:N+1)'; %求解方程
       DH0=DP(N+1);
       DP(N+1)=0;
       sum2=0;
       sum3=0;
       for i=1:N
           sum2=sum2+abs(DP(i));
           sum3=sum3+P(i);
       end
       EP=sum2/sum3; 
       P=P+0.4*DP;    %下山法，保證函數的絕對值穩定下降，加快收斂速度
   end
    H0=H0+0.3*DH0;
end
%dangliangthick;        %求HOIL

figure(1);
plot(H,'r-');    %劃出膜厚形狀曲線
hold on;
plot(P,'b-');    %劃出壓力分布曲線
hold on;
title('潤滑膜形狀和壓力分布');
axis([1 N 0 1.5]);
set(gca,'Xtick',[1:5:100],'Ytick',[0:0.2:5]);     %設置X軸的坐標從0到100，間距5;
                                                    %Y軸坐標從0到5，間距0.2

grid on;                                            %畫坐標分隔線
xlabel('Ｘ坐標點');
ylabel('潤滑膜厚/壓力值');
legend('潤滑膜厚H','壓力值P',-1);
hold off;