kill(all);
numer:true;
fpprintprec:8;retprint:false;
P_0(x):=1;
P_1(x):=1-1.0*(1-x);
P_2(x):=1.0*(-3.0*(1-x)+1.5*(1-x)^2+1);
P_3(x):=1.0*(-6.0*(1-x)-2.5*(1-x)^3+7.5*(1-x)^2+1);
P_4(x):=1.0*(-10.0*(1-x)+4.375*(1-x)^4-17.5*(1-x)^3+22.5*(1-x)^2+1);
P_5(x):=1.0*(-15.0*(1-x)-7.875*(1-x)^5+39.375*(1-x)^4-70.0*(1-x)^3+52.5*(1-x)^2+1);
P_6(x):=1.0*(-21.0*(1-x)+14.4375*(1-x)^6-86.625*(1-x)^5+196.875*(1-x)^4-210.0*(1-x)^3+105.0*(1-x)^2+1);
P_7(x):=1.0*(-28.0*(1-x)-26.8125*(1-x)^7+187.6875*(1-x)^6-519.75*(1-x)^5+721.875*(1-x)^4-525.0*(1-x)^3+189.0*(1-x)^2+1);
P_8(x):=1.0*(-36.0*(1-x)+50.273437*(1-x)^8-402.1875*(1-x)^7+1313.8125*(1-x)^6-2252.25*(1-x)^5+2165.625*(1-x)^4-1155.0*(1-x)^3+315.0*(1-x)^2+1);
P_9(x):=1.0*(-45.0*(1-x)-94.960937*(1-x)^9+854.64844*(1-x)^8-3217.5*(1-x)^7+6569.0625*(1-x)^6-7882.875*(1-x)^5+5630.625*(1-x)^4-2310.0*(1-x)^3+495.0*(1-x)^2+1);
Ns:[0,1,2,3,4,5,6,7,8,9];
Ps:[P_0, P_1, P_2, P_3, P_4, P_5, P_6, P_7, P_8, P_9];

Func(n,m,v):=block(
    index:n+1,
    h:%pi/m,buf:0,
    for i thru m-1 do buf:buf+v(i*h)*Ps[index](cos(i*h))*sin(i*h),
    return(buf*h * (2 * n + 1) / 2)
);

v(x):=x;
c:[];
for i thru 10 do(
    c:endcons(Func(i - 1, 2000, v), c)
);

R:1;
g(r,θ):=c[1]*P_0(θ)+c[2]*P_1(θ)*(r/R)+c[3]*P_2(θ)*(r/R)^2+c[4]*P_3(θ)*(r/R)^3+c[5]*P_4(θ)*(r/R)^4
+c[6]*P_5(θ)*(r/R)^5+c[7]*P_6(θ)*(r/R)^6+c[8]*P_7(θ)*(r/R)^7+c[9]*P_8(θ)*(r/R)^8+c[10]*P_9(θ)*(r/R)^9;
u (r,x):=g(r,cos(x));

u(0,0); u(0,%pi);
u(R,0); u(R,%pi);
plot3d(u(x,y),[x,0,R],[y,0,%pi],[xlabel,"радиус - r"],
    [ylabel,"угол - θ"],[legend,false],[title,"z = u(r,θ)"]);

plot2d([u(r,0),u(r,%pi/4),u(r,%pi/2)],[r,0,1],[xlabel,"радиус - r"],
    [ylabel,"u(r,θ)"],[legend,"θ=0","θ=π/4","θ=π/2"]);

plot2d([u(0.5,θ),u(0.75,θ),u(1,θ)],[θ,0,%pi],[xlabel,"угол θ"],
    [ylabel,"u(r=fix, θ)"],[legend,"r=0.5","r=0.75","r=1"]);