请问这个程序错哪里了？

我是想在Matlab中画alpha,D和Pst的三维图像，其中-1<alpha<0,0<D<1.

clc;
clear;
close;
syms x;
t=0:0.01:pi;
[D,alpha] =meshgrid(t);
beta=-0.01;%lambda=0
Q=0.50;
f=(beta/4.*D)*x^2-(2*beta./D)*x;
C=(alpha./2*D)+(3*beta/4*D)-(beta*Q/4*D.^2)+1/2;
E=(2*beta*(Q*D).^0.5)./(D.^2);
a=D*x^2+Q;
b=sqrt(D./Q)*x;
P=exp((f+E.*atan(b))).*a.^(C-0.5);
M=trapz(x,P);
N=1/M;
Pst=N*P;
mesh(D,alpha,Pst);
xlabel('D ');
ylabel(' α ');
zlabel('Pst');
axis([0 1 -1 0 0 1]);

这里是必须把x取定吗？想请大家帮我看看，应该如何修改。

看起来 最终的图片是 Pst为z轴 D，alpha 为x轴或y轴。间接参数是t 。Pst是N和P的乘积，先不看N。Pst是P，P是f,E,C的函数。f,E,C；f含D, E含D，C含alpha和D。程序运行到14行，都正确。我觉得第4行 x用syms说明它是可变的，一个x就可以画一个三维图；所以用定值比较好，title写个在x=N 时，三维图这样。
15行算不了，x需要是一列值，并且P一列一列取着算。不过这样也没有意义，这个面积看起来是要算 P这个曲面的体积值 ，Pst可能是P的单位体积归一化的一系列值，要不M改成P矩阵里的最大值吧。。。。
22行 axis范围没有曲面图啊。。。。不如自适应吧，否则应该把第5行t范围改了，但是-0.49以后公式就产生了无理数，顶多改到-0.49到pi。我还觉得 f之中含有Inf根据t取值不同不同 因为t在公式分母位置出现过，一除肯定有nan或inf啊。

clc;
clear;
close;
x=1;
%syms x;
t=-0.49:0.01:pi;
[D,alpha] =meshgrid(t);

beta=-0.01;%lambda=0
Q=0.50;
f=(beta/4*D).*x.^2-(2*beta./D).*x;
C=(alpha./2.*D)+(3*beta/4*D)-(beta*Q/4*D.^2)+1/2;
E=(2*beta*(Q*D).^0.5)./(D.^2);
a=D.*x.^2+Q;

b=sqrt(D./Q).*x;
P=exp((f+E.*atan(b))).*a.^(C-0.5);
%plot(1:1:315,P);
P(isnan(P))=1;
%M=trapz(x,P);
M = max(max(P));
N=1/M;
Pst=N*P;
mesh(D,alpha,Pst);
xlabel('D ');
ylabel(' α ');
zlabel('Pst');
%axis([0 1 -1 0 0 1]);
title(['当参数x=',num2str(x)]);