function [a,b]=Pee(r,Emax)
u=660;
e0=8.854e-12;
h=23;
j=1;
for y=0.5:0.5:60
w=atan(y/h);
syms x;
pe(1)=subs(1.5*pm(r,Emax));
pe(2)=subs(3*pm(r,Emax));
l2=subs(int(E2,x,0,u));
l1=subs(int(E1,x,0,u));
p1=roots([1/pe(1)^2+2*l2/(e0.*pe(1).*Ae(r,Emax)) 0 -1]);
p2=roots([1/pe(2)^2+2*l2/(e0.*pe(2).*Ae(r,Emax)) 0 -1]);
pn(1)=subs(int(p1(1).*E2,x,0,u))/l1;
pn(2)=subs(int(p2(1).*E2,x,0,u))/l1;
for i=3:1:10000
pe(i)=pe(i-1)+subs(pm(r,Emax)-pn(i-1)).*(pe(i-1)-pe(i-2))/(pn(i-1)-pn(i-2));
p=roots([1/pe(i)^2+2*l2./(e0.*pe(i).*Ae(r,Emax)) 0 -1]);
pn(i)=subs(int(p(1).*E2(x),x,0,u))/l1;
if abs((pn(i)-subs(pm(r,Emax)))/subs(pm(r,Emax)))<0.001
break
end
end
pend(j)=pe(i);
A(j)=sqrt(Ae(r,Emax)^2+2*pend(j)*Ae(r,Emax)*l2/e0);
b=2*pend(j)*Ae(r,Emax)*l2/e0;
E4(j)=A(j)*E(w);
p0(j)=p(1);
j=j+1;
end
for n=1:1:55
Es(n)=E4(121-n)-E4(56-n);
P(n)=p0(121-n)-p0(56-n);
end
for n=56:1:120
Es(n)=E4(121-n)-E4(n-55);
P(n)=p0(121-n)-p0(n-55);
end
for n=121:1:175
Es(n)=E4(n-120)-E4(n-55);
P(n)=p0(n-120)-p0(n-55);
end
for n=1:1:87
J(n)=-15000*P(n)*Es(n);
end
for n=88:1:175
J(n)=12000*P(n)*Es(n);
end
m=-43.5:0.5:43.5;
a=Es;
b=J;
%plot(m,Es);
%plot(m,J);
save('E720.mat','Es');
save('J720.mat','J');
end
