正方形断面のねじりのワーピング

a = 10;
b = 10;
N = 20;
x = linspace(-a,a,N);
y = linspace(-b,b,N);
[X,Y] = meshgrid(x,y);

f = @(n) (-1)^((n-1)/2)/n^3*sinh(n*pi*Y/(2*a))/cosh(n*pi*b/(2*a))*sin(n*pi*X/(2*a));
omega = 1;
Z = omega*(X*Y-32*a^2/pi^3*(f(1)+f(3)+f(5))); %n=5で打ち切り：近似
contour(X,Y,Z,'ShowText','on')
axis equal;

正解の図とは正負の分布が異なっている。