clc;clear;clf;
D=0.52;
R=D/2;
c=0.00045;
D2=0.5209;
L=1.06;
l=0.003;
u=0.096;
namda=1.2;
omiga=89;
E=100000;
v=0.49;
v0=sqrt(2*v^2/(1-v));
e=0.85;
x=0;
y=-c*e;
vx=0;
vy=0;
pq=1;
p0=0;
m=60;
n=60;
deltL=2*(L/2)/m;
deltsita=2*pi/n;
ERR=1.0e-5;
PK=zeros(n+1,m+1); aa=zeros(n+1,m+1); bb=zeros(n+1,m+1); cc=zeros(n+1,m+1);dd=zeros(n+1,m+1); ee=zeros(n+1,m+1); ff=zeros(n+1,m+1);
k=1;
while k>0
P=PK;
for i=1:n+1
for j=1:m+1
H(i,j)=1+e*(cos((i-1)deltsita))+6l*(1-v0.2)*u*omiga*R.2./(Ec.^3).P(i,j);
end
end
for i=2:n
for j=2:m