[MATLAB数学相关] 急急！哪位高手帮忙看一下我这个程序哪错了啊，有限元的

[iamwanli]

n=input('ÊäÈënµÄÖµ');
%f=inline('-(-2*x.*(1-y.^2).*(x.^2+y.^2).^(1/3).*sin(4/3*pi*(y<0)+2/3*(2*(y>=0)-1).*acos(x./(sqrt(x.^2+y.^2)+((x==0)&(y==0))))).*((x~=0)|(y~=0))+(1-x.^2).*(1-y.^2).*(2/3*(x.^2+y.^2).^(-2/3).*x.*sin(4*pi/3*(y<0)+2/3*(2*(y>=0)-1).*acos(x./sqrt(x.^2+y.^2)))-2/3*(x.^2+y.^2).^(-2/3).*y.*cos(4*pi/3*(y<0)+2/3*(2*(y>=0)-1).*acos(x./sqrt(x.^2+y.^2)))))');
%f3=inline('-(-2*y.*(1-x.^2).*(x.^2+y.^2).^(1/3).*sin(4/3*pi*(y<0)+2/3*(2*(y>=0)-1).*acos(x./(sqrt(x.^2+y.^2)+((x==0)&(y==0))))).*((x~=0)|(y~=0))+(1-x.^2).*(1-y.^2).*(2/3*(x.^2+y.^2).^(-2/3).*y.*sin(4*pi/3*(y<0)+2/3*(2*(y>=0)-1).*acos(x./sqrt(x.^2+y.^2)))+2/3*(x.^2+y.^2).^(-2/3).*x.*cos(4*pi/3*(y<0)+2/3*(2*(y>=0)-1).*acos(x./sqrt(x.^2+y.^2)))))');

%u1=inline('2*(x.^2+y.^2-2).*((x.^2+y.^2).^(1/3).*sin(4/3*pi*(y<0)+2/3*(2*(y>=0)-1).*acos(x./(sqrt(x.^2+y.^2)+((x==0)&(y==0))))).*((x~=0)|(y~=0)))-4*x.*(1-y.^2).*(2/3*(x.^2+y.^2).^(-2/3).*x.*sin(4*pi/3*(y<0)+2/3*(2*(y>=0)-1).*acos(x./sqrt(x.^2+y.^2)))-2/3*(x.^2+y.^2).^(-2/3).*y.*cos(4*pi/3*(y<0)+2/3*(2*(y>=0)-1).*acos(x./sqrt(x.^2+y.^2))))-4*y.*(1-x.^2).*(2/3*(x.^2+y.^2).^(-2/3).*y.*sin(4*pi/3*(y<0)+2/3*(2*(y>=0)-1).*acos(x./sqrt(x.^2+y.^2)))+2/3*(x.^2+y.^2).^(-2/3).*x.*cos(4*pi/3*(y<0)+2/3*(2*(y>=0)-1).*acos(x./sqrt(x.^2+y.^2))))');

%u1=inline('-2*x.*(1-y.^2).*(x.^2+y.^2).^(1/3).*sin(4/3*pi*(y<0)+2/3*(2*(y>=0)-1).*acos(x./(sqrt(x.^2+y.^2)+((x==0)&(y==0))))).*((x~=0)|(y~=0))+(1-x.^2).*(1-y.^2).*(2/3*(x.^2+y.^2).^(-2/3).*x.*sin(4*pi/3*(y<0)+2/3*(2*(y>=0)-1).*acos(x./sqrt(x.^2+y.^2)))-2/3*(x.^2+y.^2).^(-2/3).*y.*cos(4*pi/3*(y<0)+2/3*(2*(y>=0)-1).*acos(x./sqrt(x.^2+y.^2))))');
%u2=inline('-2*y.*(1-x.^2).*(x.^2+y.^2).^(1/3).*sin(4/3*pi*(y<0)+2/3*(2*(y>=0)-1).*acos(x./(sqrt(x.^2+y.^2)+((x==0)&(y==0))))).*((x~=0)|(y~=0))+(1-x.^2).*(1-y.^2).*(2/3*(x.^2+y.^2).^(-2/3).*y.*sin(4*pi/3*(y<0)+2/3*(2*(y>=0)-1).*acos(x./sqrt(x.^2+y.^2)))+2/3*(x.^2+y.^2).^(-2/3).*x.*cos(4*pi/3*(y<0)+2/3*(2*(y>=0)-1).*acos(x./sqrt(x.^2+y.^2))))');
%f=inline('-6*y-y.*(y.^2-1)+0*x');
%f3=inline('-6*x-x.*(x.^2-1)+0*y');
%u1=inline('y.*(y.^2-1)+0*x');
%u2=inline('x.*(x.^2-1)+0*y');
%g=inline('2*(x.^2+y.^2-2).*(x.^2+y.^2).^(1/3).*sin(4/3*pi*(y<0)+2/3*(2*(y>=0)-1).*acos(x./sqrt(x.^2+y.^2)))-4*x.*(1-y.^2).*(2/3*(x.^2+y.^2).^(-2/3).*x.*sin(4*pi/3*(y<0)+2/3*(2*(y>=0)-1).*acos(x./sqrt(x.^2+y.^2)))-2/3*(x.^2+y.^2).^(-2/3).*y.*cos(4*pi/3*(y<0)+2/3*(2*(y>=0)-1).*acos(x./sqrt(x.^2+y.^2))))-4*y.*(1-x.^2).*(2/3*(x.^2+y.^2).^(-2/3).*y.*sin(4*pi/3*(y<0)+2/3*(2*(y>=0)-1).*acos(x./sqrt(x.^2+y.^2)))+2/3*(x.^2+y.^2).^(-2/3).*x.*cos(4*pi/3*(y<0)+2/3*(2*(y>=0)-1).*acos(x./sqrt(x.^2+y.^2))))');
%g=inline('0*x+0*y');

f=inline('3*x.^2-6*x.*y-1-x.*y.*(y.^2-1)');
f3=inline('3*y.^2-1-6*x.*y-x.*y.*(y.^2-1)');
u1=inline('x.*y.*(y.^2-1)');
u2=inline('x.*(x.^2-1).*y');
g=inline('x.^3+y.^3-y-x');
r=zeros((2*n+1)^2,1);
%[e1,e2,C]=Merror21(n,f,f3,u1,dx,dy);h=1/n;
%[e1,C]=gbMerror21(n,f,f3,u1);h=1/n;
[C,N,T] =gbMassemble211(f,f3,g,n);h=1/n;
%[C,N,T] =bub_Massemble211(f,f3,n);h=1/n;
%e1
%e2
D=C(9*n^2+4*n+2:13*n^2+8*n+2);
c1=zeros(n^2+2*n+1,1);
c1=C(1:n^2+2*n+1);
C1=reshape(c1,n+1,n+1);

c2=zeros((2*n+1)*(n+1),1);
c2=C(n^2+n+1:3*n^2+4*n+1);
C2=reshape(c2,2*n+1,n+1);
d1=zeros(n^2+2*n+1,1);
d1=D(1:n^2+2*n+1);
D1=reshape(d1,n+1,n+1);
d2=zeros((2*n+1)*(n+1),1);
d2=D(n^2+n+1:3*n^2+4*n+1);
D2=reshape(d2,2*n+1,n+1);
p=-1:h:0;q=-1:h:0;
[P,Q]=meshgrid(p,q);C2=C2';D2=D2';
Z11=feval(u1,P,Q);
Z21=feval(u2,P,Q);
x=-1:h:1;y=0:h:1;
C1=C1';D1=D1';
[X,Y]=meshgrid(x,y);
Z12=feval(u1,X,Y);
Z22=feval(u2,X,Y);
figure(1);
mesh(P,Q,C1);hold on ;mesh(X,Y,C2);
figure(2);
mesh(P,Q,D1);hold on ;mesh(X,Y,D2);
figure(3);
mesh(P,Q,Z11);hold on ;mesh(X,Y,Z12);
figure(4);
mesh(P,Q,Z21);hold on ;mesh(X,Y,Z22);

figure(5);
x1=-1:h:1;y1=-1:h:1;
[X1,Y1]=meshgrid(x1,y1);
for k=1:n
r((k-1)*(2*n+1)+1k-1)*(2*n+1)+n+1)=C((k-1)*(n+1)+1:k*(n+1));
r(k*(2*n+1)-n+1:k*(2*n+1))=zeros(n,1);
end
r(n*(2*n+1)+12*n+1)*(2*n+1))=C(n^2+n+1:3*n^2+4*n+1);
R=reshape(r,2*n+1,2*n+1);
R=R';
contour(X1,Y1,R,40)