clc
syms  r w(r) G(r) v h r dG(r) d P E1(r) dE1(r) R

%R=1.2;% (мм) радиус пластинки
% P=35;% (мм рт ст) = p[-]-p[+]
% d=1;%(мм)радиус дырки
% v=0.45;%коэф Пуассона, сокращение в плоскости изотропии при растяжении в той же плоскости
% h=0.25;%(мм) ширина примерно пластиныы

E1(r) =5*r^2;
dE1(r)=diff(E1,r,1);
G(r)=E1(r)/(2*(1+v));
dG(r) = diff(G,r,1);


eqn = E1(r) * ( diff(w(r), r, 3) + 1/r * diff(w(r), r, 2) -...
    1/r^2 * diff(w(r), r) ) + diff(E1(r), r) *...
    (diff(w(r), r, 2) + v/r * diff(w(r), r)) == (6*P/h^3 ) *...
    (r - d/r) * (1-v^2) - diff(E1(r), r) *3*P/(5*h* dG(r) )*...
    (1 + v + d/r^2 * (1 - v ))

y2= diff(w,r,2);
y1= diff(w,r);
cond1 = [y2(d) == 6*P/(5*dG(r)*h), y1(R)==3*35/G(R)*(d/R - R),w(R)==0]
cond2 = [y2(d) == 6*P/(5*dG(r)*h), w(R)==0,-y2(R)-v*y1(R)/r+h^2/...
    (10*dG(R))*(-12*P/h^3)+6*(1+v)*v*P/(5*E1(R)*h)==0]
solw1(r) = dsolve(eqn, cond1);
solw2(r)=dsolve(eqn, cond2);


s11(r) = subs(solw(r),P,35);
s22(r)=subs(s11(r),d,1);
s33(r)=subs(s22(r),v,0.45);
s44(r)=subs(s33(r),h,0.25);

 r=d:0.001:R;
 func=s4(r);
 func2=s44(r);

 plot (r,func)
plot (r,func2)
%
% z(x,y) = (s4(sqrt(x.^2 + y.^2)))
%  for fi=0 :(pi/3):2*pi
%      for r= 0.1:0.1: 1.2
%
%        x = [x r.*cos(fi)];
%        y=[y r.*sin(fi)];
%      end
%  end
%
% [X, Y] = meshgrid(x, y);
% z =double(s4(sqrt(X.^2 + Y.^2)));
%  surf(x, y, z);
%              colormap(cool);
%              shading flat;